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I. INTRODUCTION 



The economics associated with ship operations have 
necessitated an examination of the losses associated with 
the motion of an automatically steered ship in a seaway. 

Four major areas -where fuel losses occur during the 
operation of a ship have been identified [Ref. 1, 2]. These 
areas on existing steam/diesel tankers are shown below: 

• Power plant and auxiliaries 

• Propeller efficiency 

• Hull resistance 

• Steering and navigation 

An optimized autopilot design would provide effective 
steering control with associated cost savings due to 
reducing fuel consumption. 

An appropriate computer model which represents the ship 
is necessary for studies leading to appropriate controller 
design. Chapter 2 introduces the development of two of these 
models . 

Chapter 3 addresses the formulation of a performance 
criterion which represents the added drag due to steering of 
the ship. 

Using the equations of motion as a model of the ship and 
a function minimization subroutine we proceed to the 
controller design for regular seas (deterministic model for 
the seaway) in Chapter 4, and for irregular seas (nondeter- 
ministic model) in Chapter 5. The function minimization 
subroutine used was BOXPLX and was programmed by R. R. 
Hilleary of the Naval Postgraduate School Computer Center 
[Ref. 3]. It will find the minimum of any arbitrary func- 
tion, linear or nonlinear, subject to explicit constraints 
of the variables or implicit constraints on functions of the 
variables . 
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Chapter 6 introduces another function minimization 
subroutine appropriate for onboard use. 

An adaptive control, which updates the controller param- 
eters when the environmental conditions or the ship’s course 
change, is studied in Chapter 7. 

Conclusions drawn from these experiments and recommenda- 
tions for future studies are addressed in Chapter 8. 
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II. COMPUTER MODELS OF THE SHIP 



A nontrivial part of any control problem is modelling 
the process. Thus, an appropriate computer model which 
represents the ship is necessary. The best representation of 
the ship's steering dynamics is a Taylor's series expansion 
of the force and moment relationships around a selected 
steady state operating point. The equations obtained in this 
way are known as the equations of motion [Ref. 4], and the 
formulation in the computer program is indicated in Appendix 
A. This computer program was developed by using known avail- 
able data for the SL-7 cont ainership and by implementation 
of the scheme in Figure 2.1 [Ref. 5]. 

In this scheme the function minimization subroutine is 

fed by the yaw error l y and rudder angle (5 , computes the 

r 6 

performance criterion J and adjusts the controller free 
parameters in order to minimize J. 

A second model for the ship- steering dynamics represen- 
tation is the Nomoto model. Figure 2.2 indicates the second 
and third order Nomoto transfer functions while Figure 2.3 
indicates the appropriate scheme used for obtaining these 
models from the equations of motion. Appendix A includes 
the computer program used for the Nomoto third order model 
determination . 

A yaw command is applied as input in the scheme in 
Figure 2.3 and the difference of the signals Ip^ and i- s 
fed to the function minimization subroutine which attempts 
to adjust the free parameters of 'the Nomoto plant in order 
to minimize the performance criterion J. 

Simulation runs indicate that the resulting Nomoto 
models are obtained with resulting J close to zero. 
However , in this study . the equations of motion 



13 




Figure 2.1 Controller Optimization Scheme 
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Figure 2.2 Nonioto Transfer Functions 
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Figure 2.3 Nomoto Model Determination Scheme 

representation was adopted because the system is dynamic 
and use of the Nomoto model representation implies addi- 
tional computer use. On the other hand, frequency domain 
studies were carried out using the Nomoto representation 
since this representation is easier handled. 
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III. AN ADEQUATE PERFORMANCE CRITERION 



A. CRITERION BASED ON TRUE ADDED RESISTANCE 

The performance criterion which characterizes propulsion 
losses due to steering may be shown to be that derived from 
excess power consumption per unit distance caused due to 
steering [Ref. 1, 6]. The added resistance due to steering 

can be related to the surge or thrust equation where the 
total instantaneous surge relevant to steering is 

AX=[m*( p/2)LAX; r ]y/- + l/2[(p /2)AX^]v* (3.1) 

. *1/2[(P/2)AX^ U = ] 5 2 



where 



m = 

P- 

L = 
A = 
U = 




X 



ft 



X 



vv 



mass of ship 
density of sea water 

ship's length between perpendiculars 
L 2 

ship's water speed 
sway velocity 
yaw rate of ship 
rudder angle 

= force coefficient due to yaw/ sway (posi- 
t ive ) 

= force coefficient due to rudder angle 
(negative ) 

= force coefficient due to sway 



Since the sway velocity of the ship is small we can 
neglect the term which includes the square of the sway 
velocity in the previous equation. From this the mean surge 
relevant to steering may be written as 
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AX=[m*( p/2)LAX' r ] (u g T a /2)cos( p> - ip ) 
- [(P/2)AX^U ! ]( 5 ! /2) V r 



(3.2) 



Where 



u = amplitude of 
0 

- amplitude of 

£ 

( jq - amplitude of 
^0 — ^0 - phase diff 

V r 

yaw rate 



sway ve 
yaw rat 
rudder 
erence 



locity 

e 

angle 

between 



sway 



and 



A performance criterion for added resistance due to 
steering may be formulated as 



J= lim (1/2T )/ (- a vr+7U 2 5 2 

T—*co Jn 



J o 

where Q and "V are constants 



(3.3) 



Accurate knowledge of the nonlinear coefficients and 

is required for the accuracy of such a criterion. In 
addition the criterion itself suffers from the disadvantage 
that sway velocity measurements are not available. 
Normalizing the last equation the performance criterion will 
be 



J =lim ( 1/ ZT ) / ( -X vr+ (5 2 ) dt 

norm T-+oo Jq 

where A= { 2 [m+ ( p / 2 ) LA] X ^ } / [ ( p /2)X^ 



(3.4) 



U 2 ] 



Table I indicates the 
range of speed of the ship 



values 
studied . 




for the 



operating 
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TABLE I 

Weighting factor 



Ship's speed (knots) 
16 
23 
32 




21.5350 
10.4215 
5 .3900 



B. CRITERION BASED ON APPROXIMATE ADDED RESISTANCE 

Empirical criteria based on an approximation to added 
resistance may also be derived. A semiempirical criterion 
for measuring the relative performance was developed 
[Ref. 7], based on the assumption of small amplitude oscil- 
lations around the steady-state pivot point of the ship 
during yawing at the ship/ steering system natural frequency. 
This may be extended and an alternative criterion for added 
resistance will be 



T 



J= lim 

1 —* 



(1/2T) 

CO 



5 2 ) 

0 



dt 



(3.5) 



where 



\=A(x)= [2m(l*X^)(OP/LXJ]/[( P/2)LX^U 2 ] 

\ r = [ ( p /2)LAX; r ]/m 

OP = distance from center of gravity to pivot 
center 

CO - natural frequency (closed loop ship 
steering control) 

y\ = ship's perturbation yaw angle 



The values of A as a function of ship's speed are given 
by Table II. 

A closed loop system natural frequency CO of around 0.05 
rads per sec has the potential to attenuate the effects of 
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TABLE II 

Weighting factor 

Ship’s speed (knots) 

23 
32 

seaway disturbance in the range of encounter angles where 
added resistance due to steering is important [Ref. 6]. The 
weighting factor for the operating range of the ship is 
shown in Table III. 



X 



X 



6 , 720 
3,251 
1,680 



TABLE III 
Weighting factor 



Ship ' s 



speed 

16 

23 

32 



(knots ) 



X 

X 

16 .796 
8 . 128 
4.2 



Equation 3.5 is used as a performance 
study. It is an approximation but it 
onboard use since ship's, perturbation 
rudder angle (5 are measurable. 

C. WEIGHTING FACTOR STUDY. 

The weighting factor X given by Table III used in equa- 
tion 3.5, plays an important role in terms of the optimal 
controller parameters determination. Some investigation is 
necessary in order to verify the accuracy of the results, 
since the values of X of Table III are determined based on 
the assumption that the closed loop system's natural 
frequency is around 0.05 rads per sec [Ref. 1]. Frequency 
domain techniques were used for this purpose. Using the 
Nomoto third order model representation of the ship and 
available controller parameters from Chapter 4 for sea state 



criterion for this 
is convenient for 



yaw an 



gle \p 



and 
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4, encounter frequency 1.5 rads per sec, encounter angle 
150° and ship's speed 23 knots we found that the closed loop 
bandwidth of the system is 0.04 rads per sec, as indicated 
in Figure 3.1, which is not close enough to 0.05 rads per 
sec . 

For the same sea conditions and ship speed, with the 
assumption that the closed loop natural frequency of the 
system is not 0.05 rads per sec but 0.04 rads per sec, a new 
value X=5.734 was obtained and the frequency domain tech- 
niques result in a new bandwidth 0.035 rads per sec as is 
indicated in Figure 3.2. 

Clearly, the values of X given by Table III and used in 
this study are not the best. Unfortunately, since the full 
hydrodynamic coefficients of the SL-7 cont ainership are not 
known we can't develop the surge equation and thus it is 
still impossible to determine accurate values for the 
weighting factor X. 
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( »Jlll ( til V l|*l‘ K ) 




Figure 3.1 



Closed Loop Bandwidth for 



X 



=8 . 128 
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( »•* IM ( llc'i 1 1 H‘ | v ) 




Figure 3 . 2 



Closed Loop Bandwidth for 



X 



= 5 . 734 
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IV. REGULAR SEAS 



CONTROLLER DESIGN 



We have already defined a suitable and sufficiently 
accurate ship computer model and the system's performance 
criterion, as well. The remaining task is to determine a 
representation of the external disturbances imparted to the 
ship by the sea, before the system's performance in a seaway 
can be evaluated. A correct model of the seaway itself is 
essential to representative modeling of forces and moments 
exerted on the ship by it . 

At this point we will use the regular sea model as sea 
representation. The properties of regular seas are well 
defined. The wave crests are assumed to be straight, infi- 
nitely long, parallel and equally spaced with constant wave 
height. The waves progress in a direction perpendicular to 
the crest line at a uniform velocity. However the sea is 
never regular. It is a random phenomenon where waves are 
continually changing in height, length and breadth [Ref. 8], 

The forces exerted by the regular sea have the form 



F=OJ 0 R/cos(OJ e t+T5' / . ) 



(4.1) 



wh 



The corre 
indicated in 
The excit 
cies and enco 
generator pro 



ere CJ = significant wave height 
R^* = exciting force 
(jJq- encounter frequency 
$•= phase angle 

spondence between sea state and wave 
Table IV [Ref. 9] . 

ing forces R • for different encounte 
unter angles were obtained from the 
gram [Ref . 10] . 



height is 

r frequen- 
sea state 
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TABLE IV 

Sea state vs Wave height 



Sea state 
(Beaufort scale) 
0 
1 
2 

3 

4 

5 

6 

7 

8 
9 



Range for wave height 



0 . 32 

0.65-0.98 

1.96-3.28 

3.28-4.92 

6.56-8.20 

9.84-13.1 

13.1- 18.2 

18.2- 24.6 
23.0-32.9 



Appendix B indicates the regular seastate formulation in 
the FORTRAN program used for obtaining the controller param- 
eters . 

The controller used in the entire study has one pole-one 
zero and the form is indicated in Figure 4.1. This 
controller seems to have the best performance in calm waters 
and in seaway [Ref. 5] . 




The optimized controller parameters and the cost 0 for 
23 knots speed, sea states 4-6-7-9, different encounter 
angles and various encounter frequencies are indicated in 
Tables V, VI, VII and VIII. 
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Studying the Tables V through VIII we can draw the 
following conclusions: 

• For a particular encounter angle and encounter 
frequency the higher the sea state the higher the cost. 

• For the same sea state the cost becomes smaller for 
higher encounter frequencies . 

• For encounter frequency 0.2 rads per sec the maximum 
cost occurs at 60° encounter angle for all tested sea 
states . 

• For 0.6 and 0.75 rads per sec encounter frequency the 
maximum cost occurs at 120° encounter angle for all 
tested sea states. 

• For 1.5 rads per sec encounter frequency the maximum 
cost occurs at 90° encounter angle for all tested sea 
states . 

• For 0.4 rads per sec encounter frequency the maximum 
cost occurs at 90° encounter angle for sea states 4,6 
and 7 while at sea state 9 the maximum cost occurs at 
60° encounter angle. 

Appendix C provides the computer program necessary to 
achieve the system's response. Some typical responses are 
indicated in Figures 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9. 
It is obvious that as the sea state goes to heavier seas the 
rudder and yaw perturbations become larger. 

An attempt to determine how accurate the controller 
parameters must be for a particular situation, leads to the 
conclusion that high accuracy isn't required. Keeping two 
parameters fixed each time and vary the third we can see 
(Figures 4.10, 4.11, 4.12, 4.13, 4.14) that the cost doesn’t 
change appreciably in the vicinity of the actual value. 

Figures 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9 indicate 

that the yaw and rudder excursions are less than 1° . This 
just seems strange, though it may be because of the opti- 
mization of the filter. We tried to investigate that by 
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encounter 



using the optimal filter for sea state 9, 
frequency 1.5 rads per sec, encounter angle 120° and run it 
in sea state 4, keeping the same encounter frequency and 
angle. The yaw and rudder excursions, even if they became 
larger, remained less than 1°. The parameters of those two 
filters are close and the reason might be the flatness of 
the cost surface. Second attempt led to more interesting 
'results. Using the same filter and run it in sea state 4, 
encounter frequency 0.4 rads per sec and encounter angle 
060°, the system becomes unstable (Figures 4.15, 4.16). 
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TABLE V 



Optimal Controller Parameters for Regular Sea 

Sea State 4 



Encounter Frequency 0.2 rads per sec 



encount er 
angle (degrees ) 



K1 



T1 



T2 



0 


0 . 5221067 


66.3312231 


12.8332741 


0 . 609E-33 


30 


1.0488701 


61.9309387 


15 . 9266357 


2.4819 


60 


1.2036362 


54 .5533295 


16 . 0245972 


3 .612223 


90 


1.3178062 


49 . 7426453 


14.8329315 


2 . 703524 


120 


1.3984699 


46 .9797058 


13 .9525757 


1 . 355578 


150 


1.4502153 


45.4263306 


13 .3351599 


0 . 3567196 


180 


0.7195223 


25.2119598 


14. 1219782 


0 . 345E-28 


encounter 


Encounter F 
K1 


requency^O . 4 


rads per sec 
T2 


J 


gle (degree 
0 


s) 

0.5221067 


66 . 3312231 


12 .8332741 


0 .517E-35 


30 


0 . 7234516 


74.7846533 


47.9879893 


0 . 54673 


60 


0.8730211 


92 . 8420868 


50.0454407 


1.194746 


90 


2.8910360 


28 .9679871 


10.0176315 


2 . 536161 


120 


2.6232796 


27.9702454 


8.8327761 


0.2349457 


150 


2 .6408129 


27 . 7937927 


8 .6838455 


0 . 0545772 


180 


0 . 7195223 


25.2119598 


14.1219782 


0.536E-32 


encounter 


Encounter F 
K1 


requency^O . 6 


rads per sec 
T2 


J 



angle (degrees ) 



0 


0.5221067 


66.3312231 


12 .8332741 


0 . 453E-31 


30 


2 . 0506487 


1.2107763 


5 .5998001 


0 . 0028366 


60 


1 .0504951 


0 .2329917 


19 .0010876 


0 . 0031525 


90 


2. 1496201 


1 . 2607498 


5 .3318434 


0 . 0790777 


120 


2 . 1221027 


0 . 9715905 


6 . 9064713 


0 . 0796856 


150 


1.9617786 


0.8480816 


8 .3953667 


0.0341976 


180 


0.7195223 


25.2119598 


14. 1219782 


0 . 147E-39 



Encounter 
encounter K1 

angle (degrees ) 



Frequency 0 



7 5 rads 



0 

30 

60 

90 

120 

150 

180 



0 . 

2 

2 

2 

1 

1 

0 



5221067 

4342451 

3829517 

0794163 

9938784 

9387684 

7195223 



66 .3312231 
0 . 8039263 
0.8094406 
0.4690621 
0.2545378 
0.2628353 
25 . 2119598 



per sec 
T2 



12.8332741 0 

7.4397717 0 

8.6272535 0 

16.9594727 0 

22.4892426 0 

23.8461609 0 

14.1219782 0 



7 HE- 3 5 
0017128 
0042339 
0442135 
1164415 
0379958 
312E-23 



Encounter 

encounter K1 

angle (degrees) 



Frequency^! 



5 rads 



per 

T2 



sec 



J 



0 


0 . 5221067 


66.3312231 


12 .8332741 


0 . 143E- 35 


30 


1 .8784037 


0.6553955 


5 . 4344263 


0 . 0000847 


60 


2.4894753 


0.5992758 


5 . 1160698 


0 . 0014724 


90 


2 . 5899792 


0 . 6663168 


3 . 1813316 


0 . 0241299 


120 


1 .8784037 


0 .5395500 


10.4344263 


0 . 0028635 


150 


2.4478331 


0.5559916 


5.9208412 


0.0015984 


180 


0.7195223 


25.2119598 


14. 1219782 


0 .451E-28 
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TABLE VI 



Optimal Controller Parameters for Regular Sea 

Sea State 6 



Encounter Frequency 0.2 rads per sec 



encounter 
angle (degrees) 

0 0 

30 1 

60 1 

90 1 

120 1 

150 1 

180 0 



K1 



5221067 

0966187 

2665367 

3603411 

4170542 

4533644 

7195223 



66 

62 

55 

49 

46 

45 

25 



T1 



3312231 

3766327 

1072540 

3086395 

4313202 

2830963 

2119598 



12 

19 

19 

16 

14 

13 

14 



T2 



8332741 

3941193 

2263336 

7315979 

7203903 

6198730 

1219782 



Encounter Frequency 0.4 rads per sec 
encounter K1 Tl T2 

angle (degrees ) 



0 
30 
60 
90 
120 
15 0 
180 



0 

0 

1 

2 

2 

2 

0 



5221067 

7536127 

1621914 

8681517 

6123600 

6361399 

7195223 



66 . 3312231 
12.4472111 
1.7206783 
28.5688019 
28.0108032 
27.8540497 
25.2119598 



8332741 
8151121 
6173820 
11. 3736725 
9 . 1387405 
7665482 
1219782 



12 

4 

5 



8 

14 



Encounter 

encounter K1 

angle (degrees) 



Frequency ^0 . 6 rads per^sec 



0 

8 

11 

9 

5 

1 

0 



322E-33 

622314 

84099 

270265 

004774 

396385 

678E-28 



0 

2 

4 

7 

0 

0 

0 



122E-31 

3798 

143427 

7779746 

8634676 

2140626 

345E-25 



0 


0 . 5221067 


66 . 3312231 


12 . 8332741 


0 . 806E-35 


30 


2 . 0449467 


1.2405663 


5 .5909785 


0 . 0113364 


60 


1.0594482 


0. 1432046 


18 . 5202637 


0 .0126108 


90 


2 . 1537333 


1. 1794500 


5 .9987049 


0 . 2766824 


120 


2 . 1456242 


0.9514294 


7 . 3327799 


0 .3125796 


15 0 


1.9752455 


0 .8255181 


8 .5351868 


0 . 1364259 


180 


0.7195223 


25.2119598 


14. 1219782 


0 . 112E-23 



Encounter 
encounter K1 

angle (degrees ) 



Frequency ^0 . 75 rads 



per 

T2 



sec 



0 


0.5221067 


66 . 3312231 


12 .8332741 


0. 753E-32 


30 


2.4413719 


0.7599964 


7 . 3472633 


0 . 0068484 


60 


2.4142313 


0.7818718 


8 . 5500679 


0.0169234 


90 


2.0592680 


0.3840635 


19 . 2782745 


0 . 1523162 


120 


2 . 0695038 


0.2843146 


22.5692444 


0.4640285 


15 0 


1.9513931 


1.2351496 


23.5140228 


0 . 1518951 


180 


0.7195223 


25 .2119598 


14. 1219782 


0 . 691E-28 




Encounter Frequency 1.5 r 


ads per sec 




ncount er 


K1 


Tl 


T2 


J 



angle (degrees) 



0 


0 . 5221067 


66 . 3312231 


12 . 8332741 


0 . 321E-33 


30 


1 .7606564 


0.4668925 


12 . 3475494 


0 . 0003390 


60 


2.5002985 


0.5747030 


5 . 3262844 


0 . 0058868 


90 


2 . 6015043 


0 . 6543975 


3.2408133 


0 . 0946725 


120 


2 . 1809998 


0.4698086 


10 .9160814 


0 .0114672 


150 


2.4302263 


0.5666089 


6.0247307 


0.0063947 


180 


0 .7195223 


25 .2119598 


14. 1219782 


0.344E-23 
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TABLE VII 



Optimal Controller Parameters for Regular Sea 

Sea State 7 



Encounter Frequency 0.2 rads per sec 



encounter 
angle (degrees) 


K1 


T1 


T2 


J 


0 


0.5221067 


66.3312231 


12.8332741 


0 . 911E-35 


30 


1 . 1839037 


68 .3722687 


27 .3912964 


19 . 62723 


60 


1 .3688898 


65 .5752258 


30 . 9693604 


24.53816 


90 


1.4492817 


53 . 9667511 


23 .6652069 


20 .31325 


120 


1.4652090 


46.2668457 


17 . 1322327 


12 . 58828 


150 


1.4653692 


44.8378906 


14 . 1106033 


4.036665 


180 


0.7195223 


25 . 2119598 


14 . 1219782 


0. 912E-32 


Encounter Frequency 0.4 


rads per sec 




encounter 
angle (degrees ) 


K1 


T1 


T2 


J 



0 


0.5221067 


66 . 3312231 


12 .8332741 


0 . 341E-30 


30 


0.6452138 


38.7243542 


11.6571761 


0.0524352 


60 


2.2381306 


1.0741718 


9 . 7088852 


10.22498 


90 


2 . 7780743 


30.6558838 


16.9282227 


11.95129 


120 


2 . 5894566 


28 . 1353302 


10 . 0156784 


2 . 098495 


150 


2.6306906 


27 .8760681 


8 . 9550962 


0 . 6211755 


180 


0.7195223 


25 .2119598 


14. 1219782 


0. 176E-35 



Encounter 

encounter K1 

angle (degrees) 



F requenc|_^0 



6 rads 



per sec 
T2 



J 



0 


0.5221067 


66 . 3312231 


12 .8332741 


0 . 157E-31 


30 


2 . 0501146 


1 . 2125063 


5 .5930214 


0 . 0346353 


60 


1 . 0624285 


0 . 1251755 


18.0873718 


0 . 0386276 


90 


2 . 1488247 


0 .9770966 


8 . 1855469 


0 .6328694 


120 


2 . 1970577 


0 .8763-103 


8.4206886 


0 .9178923 


150 


2.0124264 


0.8064904 


8.9983368 


0.4150548 


180 


0.7195223 


25.2119598 


14.1219782 


0 . 44 IE- 28 




Encounter F 


requency 0.75 


rads per sec 





encounter 
angle (degrees) 






T2 



0 


0.5221067 


66 . 3312231 


12 .8332741 


0.238E-28 


3.0 


2.4456072 


0 .8701959 


7 . 6628313 


0.0209509 


60 


2.4188919 


0.8008499 


8 .6116581 


0 . 0517269 


90 


1 . 9072247 


0 . 1301596 


28 .8665771 


0.4958954 


120 


2.2209492 


0.3453745 


22.8497772 


1.398130 


150 


2.0318069 


0.2577103 


23.6702576 


0.4641950 


180 


0.7195223 


25 . 2119598 


14 . 1219782 


0.542E-21 




Encounter F 


requency 1 . 5 


rads per sec 




encounter 


K1 


T1 


T2 


J 


gle (degrees ) 








0 


0.5221067 


66 .3312231 


12 .8332741 


0 .7 12E- 32 


30 


2.0834208 


0.4371362 


15 . 1762695 


0 . 0010384 


60 


2.4635468 


0 .5746776 


5 .3325920 


0 .0180051 


90 


2 .6105270 


0 . 6474890 


3.4996614 


0 . 2765892 


120 


2 . 1780138 


0.4588630 


11.4343033 


0 . 0352355 


150 


2.4314880 


0.5580992 


6.0874767 


0.0195939 


180 


0 .7195223 


25.2119598 


14. 1219782 


0.413E-26 
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TABLE VIII 



Optimal Controller Parameters for Regular Sea 

Sea State 9 



Encounter Frequency 0.2 rads per sec 



encounter 
angle (degrees) 



K1 



T1 



0 


0.5221067 


66.3312231 


30 


1 . 2344990 


97 . 2906342 


60 


1.5019569 


77 . 3820190 


90 


1.5020103 


77 . 2400513 


120 


1.5420017 


51.2350922 


150 


1.4879055 


44.4281464 


180 


0.7195223 


25.2119598 



Encounter 
encounter K1 

angle (degrees ) 



12 

53 

45 

45 

23 

15 

14 



Frequency^ . 4 rads 



T2 



,8332741 

3980713 

6973572 

6112671 

8943634 

2255249 

1219782 

per sec 
T2 



0 

31 

30 

30 

21 

8 

0 



493E-31 

88802 

87321 

87321 

75188 

599576 

810E-39 



J 



0 


0.5221067 


66 . 3312231 


12.8332741 


0. 743E-33 


30 


0.1618259 


95.8238471 


18.3397064 


8 . 103665 


60 


2.3103380 


0.9894915 


12.5459442 


20 . 30595 


90 


3 . 3482129 


2.6129665 


4.7000313 


7 . 494904 


120 


2.5573978 


28 . 7854462 


12 . 1667938 


3 . 113852 


150 


2.6185656 


27.9137726 


9 . 3553696 


1.324692 


180 


0.7195223 


25 . 2119598 


14 . 1219782 


0 . 534E-33 




Encounter Fr 


equency 0 . 6 


rads per sec 





encounter 
angle (degrees ) 



0 


0.5221067 


66.3312231 


12 .8332741 


0 .478E-31 


30 


2 . 0493793 


1.2093735 


5 . 6404839 


0 . 0820505 


60 


1.1003008 


0 . 1920759 


18 . 7369995 


0 . 0919840 


90 


2 . 0454016 


0 . 5407953 


15 . 4096680 


1 . 031631 


120 


2.2776527 


• 0.7766371 


10.4350891 


2 . 069717 


150 


2 . 0829115 


0.7725539 


9 .8721237 


0.9771649 


180 


0.7195223 


25.2119598 


14.1219782 


0 . 117E-29 



Encounter 
encounter K1 

angle (degrees ) 



Frequenc^O 



75 rads 



per 

T2 



sec 



0 


0 . 5221067 


66 . 3312231 


12 . 8332741 


0 . 872E-35 


30 


2.4433956 


0.8601446 


7 . 7009745 


0 . 0497640 


60 


2.4255047 


0.7809458 


8 . 7044067 


0. 1226490 


90 


2.2933350 


0.4751374 


16 .7149658 


0.7522082 


120 


2.3928595 


0.4103206 


22.9515076 


3 . 146294 


150 


2.1456175 


0.2891768 


23 .5924225 


1.097565 


180 


0.7195223 


25.2119598 


14. 1219782 


0. 611E-26 




Encounter Frequency 1.5 


rads per sec 




ncount er 


K1 


T1 


T2 


J 


;le (degrees ) 








0 


0.5221067 


66 .3312231 


12 . 8332741 


0 . 242E-33 


30 


2 . 0808840 ' 


0.4744592 


16 . 3423157 


0 . 0024732 


60 


2.4694090 


0.5775681 


5 .4030972 


0 . 0042757 


90 


2.6412249 


0 . 6233797 


3 . 9802380 


0 . 6127556 


120 


2 . 1955976 


0 . 4414446 


12 . 1491365 


0 . 0084519 


150 


2.4388838 


0 . 5490999 


6 . 1774960 


0 . 0046688 


180 


0 .7195223 


25.2119598 


14.1219782 


0 .420E-30 
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Figure A. 3 Rudder vs Time, Sea State A. 
Encounter Frequency 1.5 rads per sec, Encounter Angle 120 
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Figure A. A Yaw vs Time, Sea State 6. 

Encounter Frequency 1.5 rads per sec, Encounter Angle 120 
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Figure 4.5 Rudder vs Time, Sea State 6. 
Encounter Frequency 1.5 rads per sec, Encounter Angle 120 
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Figure 4.7 Rudder vs Time, Sea State 7. 
Encounter Frequency 1.5 rads per sec, Encounter Angle IzU 
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Figure A. 9 Rudder vs Time, Sea State 9. 
Encounter Frequency 1.5 rads per sec, Encounter Angle 120 
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Figure A . 10 Cost vs Kl,Sea State A. 

Encounter frequency 1.5 rads per sec, Encounter Angle 120 
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Figure A. 12 Cost vs T2,Sea State A. 

Encounter frequency 1.5 rads per sec, Encounter Angle 120 
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Figure A . 14 Cost vs T2,Sea State 6. 

Encounter Frequency 0.6 rads per sec, Encounter Angle 90 
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Figure 4.16 Rudder vs Time. Sea 4,Frequenc 
Filter for Sea State 9, Frequency 1.5 



V. IRREGULAR SEAS - CONTROLLER DESIGN 



The major characterist ic of the sea is its irregularity. 
This irregularity can be described by statistical methods by 
assuming that a large number of regular (sinusoidal) waves 
having different wavelengths,, directions, phases and ampli- 
tudes are superimposed to form the randomly varying sea. 

The presence of the irregular sea was obtained by 
coupling a sea state generator program to the FORTRAN 
program as is indicated in Appendix D. The sea state gener- 
ator program generates added mass and added inertia values 
as function of the encounter frequency and also calculates 
forces and moments imparted to the shiphull by the sea. The 
forces and moments are stored in a look up table which was 
coupled to the equations of motion. The irregular sea waves 
impinging on the ship contain the total energy density spec- 
trum composed of many frequencies and the ship responds to 
an average value of added mass and added inertia, while in 
the regular sea the added mass and added inertia were known 
for a given encounter frequency. We decided to use values 
for added mass and added inertia corresponded to encounter 
frequency 0.75 rads per sec, since the energy density is 
maximum in the vicinity of this frequency [Ref. 5]. This 
frequency gave us values representative of an average value 
for added mass and added inertia. 

The controller used for this study was the controller 
described in Chapter 4 (Figure 4.1). The optimized 
controller parameters and the cost J for sea states 4, 6, 7, 
9 and 0° ,- 30° , 60° , 90° , 120° , 150° , 180° encounter angles 
are indicated in Table IX. 

Studying the Table IX we can draw the following 
conclusions : 
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• For sea states 6, 7, 9, the higher the sea state the 

higher the cost, for every particular encounter angle. 

• Comparing costs for sea states 4 and 6 we discover some 
anomaly. The cost for a specific encounter angle in sea 
state 4 is higher than the cost for the same encounter 
angle in sea state 6. Logically, we expect higher cost 
for higher sea state. 

• The reason for this anomaly may be the method we used 
in order to obtain the added mass and added inertia 
values. The average, we consider, might not represent 
the actual average. 

Appendix E provides the computer program necessary to 
achieve the system's response. Some typical responses are 
indicated in Figures 5.1, 5.2, 5.3, 5.4. 
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TABLE IX 



Optimal Controller Paramet* 


srs for Random 


Sea 




Sea 


State 4 








encounter 


Kl 


T1 


T2 




J 


angle ( degrees 


) 










0 


0 . 6021814 


60 . 30849 


10.02579 


0 


. 342E-24 


30 


1.5121580 


89.85324 


19 . 85960 


0 


. 10719 


60 


0 . 6298036 


79 . 07199 


10 . 29221 


0 


.054196 


90 


0 . 6452737 


82 . 68692 


10 .79342 


0 


.076713 


120 


0 . 7485995 


85 .77544 


12 .37746 


0 


. 137624 


150 


0 .9101038 


92 .71379 


15 .21078 


0 


.012319 


180 


0 . 6021814 


60.30849 


10.02579 


0 


. 189E-28 




Sea 


State 6 








encounter 


Kl 


T1 


T2 




J 


angle (degrees 


) 










0 


0 . 6021814 


60 . 30849 


10 . 02579 


0 


. 172E-34 


30 


1.8743490 


61.82320 


32.22498 


0 


.044758 


60 


0 . 8662014 


90.72922 


14.44058 


0 


. 028238 


90 


0.7370305 


84 . 08502 


12 . 17457 


0 


.018541 


120 


2.6737600 


138.06650 


48 . 52447 


0 


.065028 


150 


0.4874309 


77.86977 


41.73848 


0 


.012478 


180 


0 . 6021814 


60 . 30849 


10.02579 


0 


. 767E-23 




Sea 


State 7 








encounter 


Kl 


T1 


T2 




J 


angle (degrees 


) 










0 


0 . 6021814 


60.30849 


10.02579 


0 


. 142E-34 


30 


1.7232471 


66.44597 


27 . 32489 


0 


.41237 


60 


1.8508301 


91.39410 


36 .91092 


0 


.377894 


90 


3 .7642412 


62.27244 


86 . 58169 


0 


. 258193 


120 


1.8482047 


99 . 64923 


91.28040 


0 


. 225806 


150 


0.8519831 


67.02328 


64 . 96774 


0 


.037724 


180 


0.6021814 


60.30849 


10.02579 


0 


. 811E-21 




Sea 


State 9 








encounter 


Kl 


T1 


T2 




J 


angle (degrees 


) 










0 


0 . 6021814 


60 . 30849 


10.02579 




0. 711E-35 


30 


3 . 1908160 


138 . 56101 


71.44171 




1. 306741 


60 


3.0888780 


152.01630 


72 . 66160 




0.961709 


90 


3 . 2440700 


121. 13566 


105.98480 




0 .523769 


120 


1.5461040 


111.49130 


99 . 64659 




0 .2101076 


150 


0 .3758357 


73 . 88629 


35 . 17305 




0 . 5007348 


180 


0.6021814 


60 . 30849 


10 . 02579 




0. 12 IE- 3 1 
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Figure 5.1 Yaw vs Time, Sea state 
Encounter Angle 60°. 
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Figure 5.2 Rudder vs Time. Sea State 
Encounter Angle 60°. 
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Figure 5.3 Yaw vs Time, 
Encounter Angle 




Figure 5.4 Rudder vs Time. Sea State 
Encounter Angle 60°. 



VI. MINIMIZATION SUBROUTINE FOR ONBOARD USE 



A. GENERAL 

As is mentioned earlier in Chapter 1 the function mini- 
mization subroutine used for these studies was BOXPLX. This 
subroutine will find the minimum of any function, linear or 
nonlinear, subject to explicit constraints of the variables 
or implicit constraints on functions of the variables. It 
will handle a maximum of 25 variables but can handle up to 
50 variables with user modification. 

The variables in BOXPLX are allowed to move within a 
feasible region (n-dimensional space, where n is the number 
of variables) defined by upper and lower bounds on their 
values. The choices - for upper and lower bounds for the 
parameters are based on an understanding of the function of 
each coefficient of the system. Experience indicates that 
while accurate selections of these bounds are not necessary, 
intelligent selection of these as well as the starting 
values (guesses) can considerably reduce the computer number 
of trials needed for solution convergence. This conclusion 
was drawn trying .to obtain the controller parameters for 
Tables V, VI, VII, VIII in Chapter 4. The function minimiza- 
tion subroutine, when starting the minimisation process with 
arbitrary chosen guesses, required more than 100 trials for 
convergence while by choosing guesses close to the optimal 
parameters required more than 50 and less than 100 trials. 
Considering that every trial lasted 600 seconds (10 minutes) 
and the function minimization subroutine requires about 60 
samples (trials) before telling us it had found the minimum 
this would mean 10 hours for the control to adjust itself. 
For obvious reasons, such operation is not acceptable for on 
board use. 
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For obvious reasons , such operation is not acceptable for on 
board use. 

B. ATTACKING THE PROBLEM 

We started to investigate ways to improve this. These 
efforts include: 

• Finding a more efficient function minimization 
subroutine 

• Studying the flatness of the cost surface 

• Reducing sampling time 

C. SOLVING THE PROBLEM - 

Switching to another function minimization subroutine we 
found that the new one (ZXMWD) suffered from the same disad- 
vantages . 

The experiments carried out, more than two hundred, 
indicated that the cost surface is really flat. The BOXPLX 
after a few trials started to focus on the minimum but 
before it converged, it needed more than 50 trials, even if 
the guesses were close to the optimal. The reason is the way 
BOXPLX itself tries to find the minimum of a function of NV 
variables. It converges when the cost FE remains unchanged 
for 2*NV consecutive trials with accuracy 10 s . An effort to 
modify this termination criterion in terms of the consecu- 
tive trials was successful. Table X indicates the comparison 
between modified and unmodified BOXPLX, for sea state 6, 
encounter frequency 1.5 rads per sec and encounter angle 
120°. As we can observe in Table X the cost in each case 
remains almost the same while the trials required for 
convergence are dependent on the guesses made and the termi- 
nation criterion established. 

The value of the cost is in general the summation of 
incremental contributions for each integration step and is 
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TABLE X 



First Modification in BOXPLX 



BOXPLX 


Gues ses 


Trials 


Termination 

Criterion 


Unmodified 


Arbitrary 


100 


6 


Unmodi f ied 


Optimal 


72 


6 


Modified 


Arb i t rary 


42 


3 


Modified 


Optimal 


11 


3 


Modified 


Arbitrary 


7 


2 


Modified 


Opt imal 


1 


. 2 



Cost 



0 .01146720 
0 .01146720 
0 .01146812 
0 .01146720 
0.01157926 
0 .01146720 



therefore dependent on the total time of the simulation. 
This is important in that the optimal gain coefficients 
arrived at in this manner are not optimal for steady-state 
performance but only for the time frame covered. This should 
be adequate provided the time frame selected is long rela- 
tive to the time required for the initial condition response 
to die out. This is the reason Reid has chosen time frame 
600 seconds [Ref. 1,2]. Of course this time period is large 
and we expect steady-state behaviour faster than 10 minutes. 
Simulation studies indicate that the ship, controlled by the 
controller described in Figure 4.1, reaches the steady-state 
situation in less than 100 seconds. So, we can reduce the 
600 seconds time frame to 200 seconds, safely. This is very 
important since now the modified function minimization 
subroutine BOXPLX,. uses samples of 200 seconds long instead 
of 600 seconds, converges in less than 30 minutes which is 
reasonable for on board use. Since the value of the cost is 
dependent on the time frame taken we expect reduced cost for 
200 seconds samples long, in any case. 

As we discussed earlier the BOXPLX compares the consecu- 
tive trials with accuracy 10 s . Since we do not need so big 
accuracy a second modification is necessary. We decided to 
change the existing accuracy to 10'* . 

Table XI indicates the trials required for BOXPLX 
convergence, after the second and last modification, for sea 
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TABLE XI 

Second Modification in BOXPLX 



Guesses 



Trials Termination 

criterion 



Arbitrary 17 3 
Optimal 8 3 
Arbitrary 5 2 
Optimal 1 2 



Cost 



0.003999837 
0.003995277 
0 . 004024354 
0.003995247 



state 6, encounter frequency 1.5 rads per sec and encounter 
angle 120°. By comparing Tables X and XI we can see the 
further improvement for the function minimization subroutine 
convergence. The difference in the cost is due to the 
different time frames used. The modified function minimiza- 
tion subroutine BOXPLX is indicated in Appendix F. 
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VII. ADAPTIVE CONTROL 



A. NECESSITY OF ADAPTIVITY 

The plant, the system which is supposed to be 
controlled, is normally exposed- to a time varying environ- 
ment. The ship's speed, the wave encounter angle and the sea 
state are changing drastically during the seaway. Since 
changes encountered are not completely predictable an 
optimum preprogrammed time-varying controller is not 
possible . 

If we assume-and it is apparent from the previous 
discussion in this study-that no feasible fixed parameter 
controller provides acceptable response over the entire 
performance spectrum, it is necessary that some means is 
provided for adjusting controller parameters according to 
the sea conditions and ship's operational characteristics. 
Adaptive control is thus an effort to extend basic optimum 
control concepts to these studies. 

B. CANDIDATE ADAPTIVE SCHEMES 

Since we are looking for 1% or 2% savings in fuel cost, 
we have to feed the system with precise information. Exact 
knowledge of the sea state, the wave encounter angle and the 
ship's ground speed is vital for this purpose. Currently the 
Navy is involved in a program that will provide precision 
navigation data. Garcia [Ref. 5], provides very good infor- 
mation on this subject. 

An adaptive control scheme is indicated in Figure 7.1. 
Once the adaptive part of the scheme is set up the sea 
state, encounter wave angle and ship's speed are fed to the 
filters box by the appropriate sensors. The filters box 
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Figure 7.1 Adaptive Control Scheme 
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includes predetermined optimal filter sets for different 
discrete sea states, wave angles and ship's speed. Actually, 
it is a look up table similar to those of Tables V, VI, VII 
and VIII. The output of the filter's box is a filter which 
corresponds to discrete conditions close to those fed by the 
sensors. The function minimization subroutine accepting the 
rudder angle, yaw error and the predetermined filter set as 
initial guesses tries to obtain the optimal filter for the 
exact sea state, wave encounter angle and ship's speed. At 
the same time the plant is controlled by the controller with 
the optimal predetermined set of parameters for conditions 
close to the actual. When the function minimization subrou- 
tine reaches the minimum, it supplies the controller with a 
new set of parameters which is the optimal set for the 
present conditions. 

But, what happens if either some or all the sensors 
provide new inputs to the system? A decision device placed 
after the sensors decides whether or not the change is 
appreciable. This device compares the current conditions 
with those used to obtain the controller which currently 
governs the system. If the change is higher than some 
desired percentage then a new predetermined filter is passed 
to the controller and the function minimization subroutine 
tries to find the optimal controller parameters for the new 
situation . 

From the scheme of Figure 7.1 we can eliminate the 
predetermined filters device which provides a less expensive 
system. In this case the function minimization subroutine 
will need a much longer time to determine the optimal filter 
for rapid changes in course and speed, even if we assume 
that rapid changes in sea state don’t occur. Finally, the 
function minimization subroutine might work continually "on 
line" as is indicated in Figure 7.2. In this scheme a new 
controller set is obtained in every subroutine iteration and 
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Figure 7.2 On Line Adaptive Scheme 
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some constraints for the controller parameters are necessary- 
in order to avoid operation in unstable regions. Thus, the 
filter supplied by the subroutine for controlling the plant 
must be tested for characteris t ic equation roots in the 
right half S-plane. In that case a default option must 
exist, and a special device for such purposes is necessary. 
This device will permit change in the filter parameters only 
when the new set still keeps the system in a stable situ- 
ation . 

Whichever adaptive scheme we adopt, we must provide for 
manual operation for the system. Manual operation will be 
desired in the following cases: 

• Arriving ports 

• Leaving ports 

• Restricted waters 

• Avoid collision in open seas 

• Computer down 

For the first two situations, since we usually expect no 
heavy seas, the optimal filter for sea state 1 seems to be 
more suitable. Also, this filter can serve as initial condi- 
tion for the adaptive schemes described before, when we are 
leaving ports. For the rest of the situations the last 
operating optimal filter is the most appropriate. 
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VIII. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 

The principal conclusions from this study of the SL-7 
containership as they related to steering may be stated as 
follows : 

• It is evident that a control system which provides the 
ship heading and simultaneously reduces the propulsive 
losses does exist and therefore such a controller saves 
fuel. We can't conclude how much the savings are, 
since there is no reference for comparison between' the 
conventional autopilot and the autopilot which, in 
addition, holds the potential for reducing the propul- 
sive losses. The literature says that savings 1% or 2 % 
is possible. 

• An adaptive controller, that minimizes propulsion 
losses as ship characteris t ics and .environmental condi- 
tions change, may be designed using a self- opt imalizing 
technique employing a suitable performance criterion. 

• Studying every particular situation we conclude that 
the cost surface is flat and therefore accurate deter- 
mination of the controller parameters is not required. 
So, if we decide to use the adaptive control scheme of 
Figure 7.1 we don't have to store every particular 
filter in the look up table since one filter may be 
suitable for different ship characteristics and sea 
condi t ions . 

• The weighting factor X used in the performance 
criterion equation 3.5 plays an exceptional and impor- 
tant role in the optimal controller parameters determi- 
nation. However, it is obtained from studies based on 
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many assumptions and therefore isn't predicted accu- 
rately. As a consequence the controller found does not 
minimize added resistance unless the proper value of 
the weighting factor X has been used. 

• The method used to obtain the average for the added 
mass and added inertia for the irregular seas studies 
might not represent the actual average. 

• The function minimization subroutine BOXPLX after two 
modifications seems to be pretty suitable for on board 
use working as a main part of the adaptive scheme. 
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RECOMMENDATIONS FOR FUTURE STUDIES 



The following recommendations for future work to gain a 
fuller and deeper understanding of the problem are made as 
follows : 

• Some studies are necessary in order to investigate why 
part of the obtained controllers in Tables V, VI, VII, 
VIII are lag and part of them are lead filters. 

• As we stated in chapter 4 the yaw and rudder excursions 
in Figures 4.2 through 4.9 are less than 1 J . It is 
necessary to investigate the reason for that. It might 
be because of the optimization of the filter or the 
forces and moments have not been sealed properly. 

• It is necessary to find out the appropriate average 
values for the added mass and added inertia before we 
attempt further studies in irregular seas . 

• The full hydrodynamic coefficients for the SL-7 are 
necessary in order to develop and include the surge 
equation in our ship model. So far we ignore the surge 
equation and we assumed constant ship speed while in 
reality the added resistance due to steering must 
reduce the ship speed. 



63 



• By developing the surge equation in our ship model we 
may be able to determine a good value for the weighting 
factor to be used in the performance criterion. 

• Also, with the surge equation available in our model we 
should be able to calculate actual energy losses and 
savings in fuel. 

• As we stated in chapter 6 the time frames (samples) 
used by the BOXPLX must be long relative to the time 
required for the initial condition response to die out. 
Reid [Ref. 1,2] has chosen time frame 600 seconds long. 
This time frame is long and it is necessary to find out 
the sufficient time frame since the time required by 
the function minimization subroutine for convergence is 
directly proportional to the sample duration. 
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APPENDIX A 



NOMOTO THIRD ORDER MODEL DETERMINATION 

//NOMOTO JOB (XXXX,XXXX) , ’RESEARCH' , CLASS =J 

//-MAIN ORG=NPGVMl.XXXXP 

// EXEC FORTXCG , PARM . FORT= ’ OPT ( 2 ) ’ , IMSL=DP , REGION= 1024K 

//FORT. SYS IN DD * 

C IN ORDER TO PERFORM SIMULATION ONLY WHEN GAINS HAVE BEEN 

C OBTAINED CHANGE XS(-) TO X(*) AND DELETE XU ( - ) , AND XL ( ) 

DIMENSION XS(4) ,XU(4) ,XL(4) 

XS(1)=0. 1 
XS ( 2 ) = 15 . 13 
XS (3 )=15 .675 
XS (4)=9 . 014 

C XS(I) IS THE STARTING GUESS 

C XL(I) IS THE LOWER LIMIT FOR THE I ’ TH VARIABLE 

C • XU(I) IS THE UPPER LIMIT FOR THE I ’ TH VARIABLE 

XL ( 1 ) = .01 
XU ( 1 ) = 1 . 0 
XL ( 2 ) = 1 . 0 
XU ( 2 ) - 2 0 . 

XL ( 3 ) = 1 . 0 
XU (3 ) = 100 . 

XL ( 4 ) = 1 . 0 
XU (4 ) = 100 . 0 

C A DESCRIPTION OF THE FOLLOWING PARAMETERS 

C IS DISCUSSED IN BOXPLX 
R= 9 . / 13 . 

NTA= 1000 
NPR=0 
NAV= 0 
NV=4 
IP = 0 
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C THE FOLLOWING STATEMENT MUST BE CHANGED TO 
C CALL PLANT (XX) 

C IF ONLY SIMULATION IS WANTED 

CALL BOXPLX ( NV , NAV , NPR , NTA , R , XS , IP , XU , XL , YMN , IER) 
WRITE (6,25) 

25 FORMAT (IX,’ OPTIMAL GAINS ',/ ) 

DO 30 1=1,4 

30 WRITE (6 ,40)I,XS(I) 

40 FORMAT (IX, 'X( ’ ,12 , ' )= ' ,F14. 7) 

WRITE (6,77) TDIFF 
77 FORMAT (IX, ’ COST= ' ,F14. 7) 

STOP 

END 

C SUBROUTINE PLANT (XX) SIMULATES THE SHIP 
- SUBROUTINE PLANT (XX) 

COMMON TDIFF 

REAL *8 L , L2 , L3 , L4 , L5 , L6 , RXR , RXI , RYR , RYI , MZR , MZI , RX , RY 
REAL* 8 X , XDOT , Y , YDOT , U , UDOT , V , VDOT , YAW , R , RDOT , TX , TY 
REAL* 8 TIME , ETIME , XUDOT , XU , XUU , XVR , XVV , XDD , WA , WE , RZ 
REAL* 8 YV , YR , YD , YVVR , YVRR , YVVV , YRRR , YDDD , YVDOT , TZ 
REAL* 8 NV , NR , ND , NVVR , NVRR , NVVV , NRRR , NDDD , NRDOT 
REAL *8 RHO , IZ , FX , FY , MZ , XP , MASS , DELT 
REAL* 8 YAWE , YAWC , D2 , D 

REAL* 8 K , TP 1 , TP2 , Z , XI ,X2 , X3 , X4 , X5 , YAW2 , S 
DIMENSION XX (4) 

C INITIAL CONDITIONS FOR INTEGRATION 
C SIMULATION END TIME IN SECONDS 
ETIME = 600 . 

TIME = 0 . 0 
ICOUNT= 1 

C INITIALIZE THE COST FUNCTION 
TDIFF=0 . 0 

C GAIN COEFFICIENTS TO BE OPTIMIZED 
K=XX( 1 ) 

Z=XX ( 2 ) 
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TP 1=XX ( 3 ) 

TP2=XX ( 4 ) 

C X , XDOT , Y , YDOT ARE FIX COORDINATES ON EARTH 
X = 0. 0 
Y=0 . 0 
XDOT=0 . 0 
YDOT-O . 0 

C U , UDOT , V , VDOT ARE FIX COORDINATES ON SHIP 
V=0. 0 
UDOT = 0 . 0 
VDOT= 0 . 0 
YAW =0.0 
R=0 .0 
RDOT=0 .0 

C ORDERED SPEED IN FEET/ SEC 

C 38.81 FT/ SEC=23 KNOTS 
UC=38 .81 

C AT STEADY STATE ACTUAL SPEED (U) = COMMAND SPEED (UC) 
U = UC 

C D = RUDDER ANGLE 
D= 0 . 0 
L=880 . 5 
L2=L**2 
L3=L*L*L 
L4=L*L3 
L5=L*L4 
L6=L*L5 

C SEA DISTURBANCE 

C FORCES IN X,Y DIRECTION COMPUTED IN FORCES 

C MOMENTS IN Z ' 

FX=0 . 0 
FY=0 . 0 
MZ = 0 . 0 

C RXR=- . 15744D+05 

C RXI=- . 19950D+06 
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c 


RYR=0 . 52365D+04 




c 


RYI=0 . 18699D+06 




c 


MZR= - . 29870D+08 




c 


MZI= - . 3775 1D+07 






RXR= - . 50230D+04 






RXI = 0 . 127 12D+ 05 






RYR= 0 . 35290D+ 04 






RYI= - . 31909D+05 






MZR = 0 . 38826D+07 






MZI= - . 643 13D+07 




c 


RXR= 0 . 28540D+04 




c 


RXI= - . 99574D+04 




c 


RYR= - . 85441D+04 




c 


RYI=0 .39595D+05 




c 


MZR= - . 130 14D+ 08 




c 


MZI=0 . 11348D+08 




c 


RXR= - . 7 5 642D+ 04 




c 


RXI=0 . 83497D+04 




c 


RYR= 0 . 23379D+05 




c 


RYI= - . 8150 2D +05 




c 


MZR=0 . 28622D+07 




c 


MZI = - . 19388D+ 08 




c 


RXR= - .37916D+04 




c 


RXI=0 . 1638 1D+04 




c 


RYR= - . 76647D+05 




c 


RYI= - . 37685D+05 




c 


MZR- - .839 15D+07 




c 


MZI= - .53176D+07 






RX=DSQRT ( RXR** 2 + RX I * 


*2) 




RY-DSQRT ( RYR** 2 + RYI * 


*2 ) 




RZ = D S QRT ( MZR** 2 + MZ I * 
TX= DATAN2 ( RXI , RXR ) 

T Y = D ATAN 2 ( RYI , RYR ) 
TZ=DATAN2 (MZI ,MZR) 


*2 ) 



C SIGNIFICANT WAVE HEIGHT; SEA STATE 1-0.32,2-0.75,3-2.5, 
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C 4-5.0,5-7.0,6-10.0,7-17.5,8-20.5,9-27.0 
WA= 10 . 0 

C ENCOUNTER FREQUENCY .1,. 2, .3, .4, .5, .6, .75, 1.0, 1.5, 2. 5 
WE=0 . 4 

C HYDRODYNAMIC COEFFICIENTS ARE INSERTED AS PARAMETERS 
RHO =1.9876 

MAS S = ( . 0 0 44 ) * ( . 5 *RHO*L3 ) 

IZ= (0 . 00028 )*( .5 *RHO*L5 ) 

YAWE = 0 . 0 
X1=0.0 
X2 = 0 . 0 
X3 = 0 . 0 
X4 = 0 . 0 
X5 = 0 . 0 
YAW2=0 .0 
200 CONTINUE 

S = D S QRT ( U** 2 + V**2) 

C INPUT YAW COMMAND 
YAWC = 0 . 0 

IF (TIME. GE. 0.0) YAWC= ( 1 . 0/ 57 . 296 ) 

C ERROR SIGNAL TO DRIVE RUDDER (YAW ACTUAL - YAW COMMAND) 
C FOR EQUATIONS OF MOTION. 

YAWE=YAW - YAWC 
D=YAWE 
C 

C NOMOTO 3RD ORDER PLANT 
C 

C ERROR SIGNAL TO DRIVE RUDDER (YAW COMMAND - YAW ACTUAL) 
C FOR NOMOTO MODEL. 

D2 = YAWC - YAW2 
X1=(D2-X2)/TP1 
X3=K*(Z*X1+X2 ) 

X4= (X3-X5 )/TP2 

C AXIAL FORCE HYDRODYNAMIC COEFFICIENTS (SURGE) 

XUDOT =(-.0001) *(.5 *RHO*L3 ) 
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XUU = ( - 0 . 0 0 0 3 ) * ( . 5 *RH0*L2 ) 

XU= ( - 0 . 0 25 3 ) * ( 0 . 5 *RH0*L2*S ) 

XVR= ( 0 . 00 3 9 ) * ( . 5 *RHO*L3 ) 

XVV= ( - . 0012 ) * ( . 5*RHO*L2 ) 

XDD= ( - 0 . 0005 )* ( . 5*RHO*L2*S**2 ) 

C LATERAL FORCE HYDRODYNAMIC COEFFICIENTS (SWAY) 

C YV= ( - 0 . 00 7 5 8 ) * ( . 5*RHO*L2*S ) 

YR= (0 . 0023 )*(’. 5*RHO*L3*S ) 

YD= (0 . 00145 )*( . 5*RHO*L2*S**2) 

YVVR = ( 0 . 0 1 ) * ( . 5 *RHO*L3 / S ) 

YVRR= ( - 0 . 008 ) * ( . 5 *RHO*L4 / S ) 

YVVV= ( - 0 . 0 3 ) * ( . 5 *RHO*L2 / S ) 

YRRR= ( 0 . 003 ) * ( . 5 *RHO*L5 / S ) 

YDDD- (-0 . 0005 )*( . 5*RHO*L2*S**2 ) 

YVDOT = - 0 . 30908D+ 07 
YV= -0. 81271D+04 
C YVDOT= - . 36185D+07 

C YV=- . 24757D+06 

C YVDOT= - . 32890D+07 

C YV= - . 11775D+07 

C YVD0T= - .23038D+07 

C YV=- . 18267D+07 

C YVDOT=- .59800D+06 

C YV= - . 13260D+07 

C MOMENT ABOUT Z-AXIS HYDRODYNAMIC COEFFICIENTS (YAW) 

NV= ( - 0 . 00213 )* ( . 5*RHO*L3*S ) 

C NR= ( - 0 . 00 10 5 ) * ( . 5 *RH0*L4*S ) 

ND= ( - 0 . 000 7 ) * ( . 5 *RHO*L3*S**2 ) 

NVVR= ( - 0 . 0 15 ) * ( . 5 *RH0*L4 / S ) 

NVRR= ( - 0 . 008 ) * ( . 5*RHO*L5 / S ) 

NVVV= ( 0 . 0 1 ) * ( . 5 *RHO*L3 / S ) 

NRRR= (-0.006 )*( .5*RHO*L6/S) 

NDDD= ( 0 . 0 0 0 1 ) * ( . 5 *RHO *L3 * S ** 2 ) 

C NRDOT IS THE ADDED INERTIA TERM WHICH MUST BE CHANGED 
C FOR DIFFERENT ENCOUNTER ANGLE , SPEED , ENCOUNTER FREQUENCY 
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c 

c 

c 

c 

c 



c 

c 

c 

c 

c 

c 

c 



c 

c 

c 

c 

c 

c 



c 



NRDOT= ( - 0 . 00027 )* ( . 5*RHO*L5 ) 

SPEED r 23 KNOTS , ENCOUNTER ANGLE = 60, ENCOUNTER FREQ=0.75 
NRDOT= - . 2625 1D+ 12 
NR= - . 53637D+09 
NRDOT= - . 20 125D+ 12 
NR= - .-949 7 0D+ 10 
NRDOT= - . 1867 1D+ 12 
NR=- . 46860D+ 11 
NRDOT= - . 14518D+12 
NR=- . 87538D+11 
NRDOT= - . 37261D+11 
NR= - . 69856D+11 
SETS SEA STATE TO ZERO 
FX= 0 . 

FY= 0 . 

MZ = 0 . 

F X = WA*RX "D COS (WE*TIME + TX) 

F Y= WA*RY*DCO S ( WE* TIME + TY ) 

M Z = W A * R Z * D C O S (WE*TIME + TZ ) 

U ACTUAL SPEED 
UC COMMANDED SPEED 
XP = PROPELLER THRUST 
XP=-XUU*UC**2 



EQUATIONS OF MOTION 

UDOT= ( (XVR + MASS ) *V*R + XUU*U**2 + XVV*V**2 
1 + XDD*D*D + FX + XP ) / (MASS-XUDOT) 



VDOT= ( YV*V + ( YR - MA S S *U ) *R 

1 + YVRR*V*R**2 + YVW*V*V* 

2 + YRRR*R**3 + YDDD*D**3 + 
RDOT=(NV*V + NR*R + ND*D + 

1 +NVW*V**3 + NRRR*R**3 + 
WHEN TO PRINTOUT 

IF (ICOUNT.EQ.il) GO TO 50 
GO TO 300 



+ YD'-'D + YWR*V**2*R 
*3 

FY )/ (MASS-YVDOT) 
NWR*V**2*R + NVRR*V*R**2 
NDDD*D**3 + MZ ) / ( I Z - NRDOT ) 
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C CONVERT RADIANS TO DEGREES 
50 YAWDEG= YAW*57.296 
RDEG=R*57 . 296 
RDD E G = RDOT* 5 7 .29 6 
DDEG=D*57 . 296 
YAWC = YAWC* 5 7 .29 6 
ICOUNT= 1 

C TEST IF WANT TO STOP 
300 IF ( TIME . GE . ETIME ) GO TO 400 
C INTEGRATION STEP SIZE DELT 
DELT =1.0 
C INTEGRATION 

X2=X2+X1*DELT 
X5=X5 +X4*DELT 
YAW2=YAW2+X5*DELT 
U=U+UDOT*DELT 
V = V + VDOT * DELT 
R=R+RDOT*DELT 
YAW= YAW+R' V DELT 

C CONVERT SHIP TO FIXED COORDINATES ON EARTH 
XDOT = U*DCOS (YAW ) - V*D-SIN (YAW ) 
YDOT=U*DSIN(YAW)+V*DCOS(YAW) 

X=X+XDOT' v DELT 
Y= Y + YDOT' v DELT 
TIME=TIME+DELT 
ICOUNT= ICOUNT+ I 
C COST FUNCTION 

TDIFF=TDIFF+ (YAW-YAW2)** 2 
GO TO 200 
400 CONTINUE 

C WRITE (6 , 500 ) TDIFF , K , Z , TP 1 , TP2 

500 FORMAT ( ’ ’ , IX , ’ COST = ’ , F12 . 7 , 2X , ’ K =’,F10.7, 

I ’ Z = ’ , F15 . 7 , ’ TP I = ’ ,F15 .7 ,2X, ’ TP2 =',F15.7) 
RETURN 
END 



72 



APPENDIX B 



REGULAR SEASTATE FORMULATION 

//REGUGAINS JOB (XXXX , XXXX ),’ RESEARCH CLAS S = C 

//-'MAIN ORG = NPGVMl.XXXXP 

// EXEC FORTXCG , PARM . FORT= ’ OPT ( 2 ) ’ , IMSL=DP , REGION= 1024K 

//FORT. SYS IN DD - 

C IN ORDER TO PERFORM SIMULATION ONLY WHEN GAINS HAVE BEEN 

C OBTAINED CHANGE XS(-) TO X(*) AND DELETE XU ( - ) , AND XL(-) 
DIMENSION XS(3) , XU ( 3 ) ,XL(3) 

XS(1)=0 .9650610 
XS (2) = 0 .4500911 
XS(3)=5 .6194260 

C XS(I) IS THE STARTING GUESS 

C XL(I) IS THE LOWER LIMIT FOR THE I ’ TH VARIABLE 

C XU(I) IS THE UPPER LIMIT FOR THE I’TH VARIABLE 
. XL ( 1 ) = 0 . 1 
XU ( 1 ) = 4 . 0 
XL ( 2 ) = 0 . 1 
XU(2)=15 .0 
XL ( 3 ) = 1 . 0 
XU ( 3 ) = 2 5 . 0 

C A DESCRIPTION OF THE FOLLOWING PARAMETERS 

C IS DISCUSSED IN BOXPLX 
R= 9 . / 13 . 

NTA= 1000 
NPR= 100 
NAV= 0 
NV= 3 
IP = 0 

C THE FOLLOWING STATEMENT MUST BE CHANGED TO 

C CALL PLANT ( X ) 

C IF ONLY SIMULATION IS WANTED 
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CALL BOXPLX ( NV , NAV , NPR , NTA , R , XS , IP , XU , XL , YMN , IER ) 
WRITE (6,25) 

25 FORMAT (IX,' OPTIMAL GAINS ',/ ) 

DO 30 1=1,3 

30 WRITE(6 ,40)1 ,XS(I) 

40 FORMAT (IX, 'X(* ,12, ' )=' , F14 . 7 ) 

STOP 

END 

SUBROUTINE PLANT (XX) 

C SUBROUTINE PLANT (XX) SIMULATES THE SHIP 
COMMON TDIFF 
REAL" 8 L,L2 ,L3 ,L4 ,L5 ,L6 

REAL "8 X , XDOT , Y , YDOT , U , UDOT , V , VDOT , YAW , R , RDOT 
REAL* 8 TIME , ETIME , XUDOT , XUU , XVR , XW , XDD 
REAL "8 YV , YR , YD , YWR , YVRR , YVW , YRRR , YDDD , YVDOT 
REAL "8 NV , NR , ND , NVVR , NVRR , NVVV , NRRR , NDDD , NRDOT 
REAL* 8 RHO , IZ , FX , FY , MZ , XP , MAS S , DELT , MZI , WA , WE 
REAL* 8 DYAWE , YAWE , YAWC , ISE , I SR , LAMDA , D , RYR , RYI , MZR 
REAL* 8 K1 , T1 , T2 , D , X2 , DX2 , S , RX , RY , RZ , TX , TY , TZ , RXR , RXI 
DIMENSION XX ( 3 ) 

C 

C CLOSE LOOP ANALYSIS WITH FILTER 

C 

C INITIAL CONDITIONS FOR INTEGRATION 

C SIMULATION END TIME IN SECONDS 
ETIME= 6 00 . 0 
TIME=0 .0 
ICOUNT = 1 

C INITIALIZE THE COST FUNCTION 
ISE=0. 0 
ISR=0 . 0 
TDIFF = 0.0 
LAMDA=8 . 128 

C- GAIN COEFFICIENTS TO BE OPTIMIZED 
K1=XX ( 1 ) 
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T1=XX ( 2 ) 

T2=XX(3 ) 

C WRITE(6 , 1010) K1 , T1 , T2 

C1010 FORMAT (IX, 'K1 =’,F15.7,’ T1 = ’ , F 15 . 7 , ’ T2 =\F15.7) 

C X , XDOT , Y , YDOT ARE FIX COORDINATES ON EARTH 
X=0.0 
Y=0 . 0 
XDOT=0 . 0 
YDOT = 0 . 0 

C U , UDOT , V , VDOT ARE FIX COORDINATES ON SHIP 
V-0 . 0 
UDOT = 0 . 0 
VD0T=0 . 0 
YAW=0 . 0 
R=0 .0 
RDOT= 0 . 0 

C ORDERED SPEED IN FEET/ SEC 

C 38.81 FT/ SEC = 23 KNOTS 
UC= 38 .81 

C AT STEADY STATE ACTUAL SPEED (U) = COMMAND SPEED (UC) 

U = UC 

C D = RUDDER ANGLE 
D= 0 . 0 
L= 880 . 5 
L2=L**2 
L3 = L*L‘ V L 
L4=L*L3 
L5 = L~'L4 
L6=L*L5 

C SEA DISTURBANCE 

C FORCES IN X,Y DIRECTION COMPUTED IN FORCES 

C MOMENTS IN Z 
FX=0 . 

FY= 0 . 

MZ = 0 . 
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C RXR= -0 . 91037D+03 

C RXI=0 . 50869D+05 

C RYR= -0 . 20256D+04 

C RYI=0 . 18077D+06 

C MZR=- . 14310D+08 

C MZI= - . 16903D+07 

C RXR= -0 . 99047D+04 

C RXI= . 15994D+06 

C RYR= - . 64455D+05 

C RYI=0 . 61873D+06 

C MZR= . 120180+08 

C MZI=- .49204D+07 

C RXR=-0 . 32876D+05 

C RXI=0 . 25844D+06 

C RYR=- .27053D+06 

C RYI= . 90191D+06 

C MZR=0 . 11964D+09 

C MZI=0 . 24103D+08 

C RXR= - . 54639D+05 

C RXI= . 28236D+06 

C RYR= - . 28668D+06 

C RYI=0 .79670D+06 

C MZR=0 . 19925D+09 

C MZI=0 .77746D+08 

RXR= 0 . 27268D+05 
RXI= - . 71601D+05 
RYR=0 .14077D+05 
RYI= - . 28679D+06 
MZR= - . 30892D+08 
MZI= - . 53246D+08 
RX = D S QRT ( RXR' V ' V 2 + RX I * * 2 ) 
R Y = D S QRT ( R YR** 2 + R Y I ** 2 ) 
RZ = DSQRT (MZR**2 + MZI**2 ) 
TX=DATAN’2 (RXI,RXR) 

T Y = DATAN 2 ( RYI , RYR ) 



76 



TZ=DATAN2 (MZI , MZR ) 

C SIGNIFICANT WAVE HEIGHT; SEA STATE 1-0.32,2-0.75,3-2.5, 
C 4-5.0,5-7.0,6-10.0,7-17.5,8-20.5,9-27.0 
WA= 10 . 0 

C ENCOUNTER FREQUENCY .1,. 2, .3, .4, .5, .6, .75, 1.0, 1.5, 2. 5 
WE = 1 . 5 

C HYDRODYNAMIC COEFFICIENTS ARE INSERTED AS PARAMETERS 
RHO = 1 . 98 7 6 

MASS = ( . 0044 ) * ( . 5 “RHO' v L3 ) 

•IZ= (0 . 00028 )*( . 5*RHO*L5) 

YAWE= 0 . 0 
X2 = 0 . 0 
DX2 = 0 .0 
200 CONTINUE 

S=DSQRT (U**2+V**2 ) 

C INPUT YAW COMMAND 
YAWC = 0 . 0 

IF (TIME. GE. 0.0) YAWC=0.0 

C. ERROR SIGNAL TO DRIVE RUDDER (YAW ACTUAL - YAW ORDERED) 

C ( COMPENSATOR FILTER ) 

YAWE= YAW - YAWC 
DX2= (YAWE-X2 )/T2 
D=K1*(T1*DX2+X2) 

C AXIAL FORCE HYDRODYNAMIC COEFFICIENTS (SURGE) 

C XUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED FOR 
C DIFFERENT ENCOUNTER ANGLE , SPEED , ENCOUNTER FREQUENCY 
C 

XUDOT = (-.0001) *(.5 -’'RHO *L 3 ) 

XU= ( - 0 . 025 3 ) * ( . 5*RHO*L2*S ) 

XUU= ( - 0 . 00 0 3 ) * ( . 5 “'RHO“L2 ) 

XVR= ( 0 . 0 0 3 9 ) * ( . 5 “RHO"'L3 ) 

XVV= ( - . 00 12 ) “ ( . 5*RHO*L2 ) 

XDD= (-0.0005) *(.5 “RH0“'L2"'S“"2) 

C LATERAL FORCE HYDRODYNAMIC COEFFICIENTS (SWAY) 

C YV= ( - 0 . 00758 ) * ( . 5 “RHO“L2 “S ) 
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YR= (0 . 0023) * ( . 5*RHO*L3*S ) 

YD= (0 . 00145) * ( . 5*RHO*L2*S**2 ) 

YVVR= ( 0 . 0 1 ) * ( . 5 *RHO*L3 / S ) 

YVRR= ( - 0 . 00 8 )' ; '(. 5 *RH0*L4 / S ) 

YVVV= ( - 0 . 0 3 ) * ( . 5 *RHO*L2 / S ) 

YRRR= ( 0 . 003 ) * ( . 5*RHO*L5 / S ) 

YDDD= ( - 0 . 0005 ) * ( . 5*RHO*L2*S**2 ) 

C YUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED FOR 
C DIFFERENT ENCOUNTER ANGLE , SPEED , ENCOUNTER FREQUENCY 
C 

C YVDOT= (-0 . 0039 )*( . 5*RHO*L3) 

C SPEED=23 KNOTS , ENCOUNTER ANGLE= 60 , ENCOUNTER FREQ=0.75 
C YVDOT= - . 30908D+07 

C YV= - . 8 127 1D+04 

C YVDOT= -.36185D+07 

C . YV= - . 2475 7D+ 06 
C YVDOT=- . 32890D+07 

C YV= - . 11775D+07 

C YVDOT=- . 23038D+07 

C YV= - . 18267D+07 

YVDOT= - . 59800D+06 
YV= - . 13260D+07 

C MOMENT ABOUT Z-AXIS HYDRODYNAMIC COEFFICIENTS (YAW) 

NV= (-0 . 00213) *( . 5*RHO*L3*S) 

C NR= ( - 0 . 00 105 ) * ( . 5*RHO*L4*S ) 

ND = ( - 0 . 0 0 0 7 ) * ( . 5 *RHO *L3 * S ** 2 ) 

NWR : ( - 0 . 0 15 ) * ( . 5 *RH0*L4 / S ) 

NVRR= (-0.008)* (.5 *RH0*L5 / S ) 

NVVV= ( 0 . 01 ) * ( . 5 *RHO*L3 / S ) 

NRRR= ( - 0 . 0 0 6 ) * ( . 5 *RHO*L6 / S ) 

NDDD= (0 . 0001 )*( . 5*RHO*L3*S**2 ) 

C NRDOT IS THE ADDED INERTIA TERM WHICH MUST BE CHANGED 
C FOR DIFFERENT ENCOUNTER ANGLE , SPEED , ENCOUNTER FREQUENCY 
C 

C NRDOT = ( - 0 . 0 0 0 2 7 ) * ( . 5 *RHO*L5 ) 
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c 

c 

c 

c 

c 

c 

c 

c 

c 



c 



c 

c 

c 

c 

c 

c 



c 



c 



SPEED= 23 KNOTS , ENCOUNTER ANGLE = 60 , ENCOUNTER FREQ=0.75 
NRDOT= - . 2625 1D+ 12 
NR= - . 53637D-09 
NRDOT= - . 20 125D+ 12 
NR=- . 94970D+10 
NRDOT= - . 18 6 7 1D+ 12 
NR=- . 46860D+ 11 
NRDOT= - . 145 18D+ 12 
NR= - . 87538D+11 
NRDOT= - . 3726 1D + 11 
NR= - . 69856D+11 
REGULAR WAVE SEA STATE 

FX = WA*RX*DCO S ( WE "TIME + TX ) 

F Y = WA*RY*D C 0 S ( WE "TIME + TY ) 

MZ =• WA " RZ " D C 0 S ( WE "TIME + TZ ) 

U ACTUAL SPEED 
UC COMMANDED SPEED 
XP = PROPELLER THRUST 
XP=-XUU"UC""2 
EQUATIONS OF MOTION 

UDOT= ( (XVR + MASS ) "V"R + XUU*U**? + XVV"V""2 

1 + XDD"'D*D + FX + XP )/ (MASS-XUDOT) 

VDOT=(YV"'V + (YR-MASS"U)"'R 

1 + YVRR"V"'R""2 + YVW-'V""'3 

2 + YRRR "'R " "' 3 + YDDD*D**3 + 

RDOT= (NV*V + NR"R + ND"D + 

1 + NWV"'V""3 + NRRR*R""3 + 

WHEN TO PRINTOUT 

IF (ICOUNT.EQ.il) GO TO 50 
GO TO 300 

CONVERT RADIANS TO DEGREES 
50 YAWDEG= YAW"'5 7.29 6 



+ YD"D + YWPV"V""2"R 

FY )/ (MASS-YVDOT) 

NVVR " V " " 2 "'R + NVRR " V *R " 2 
NDDD"'D" "3 + MZ ) / ( IZ - NRDOT ) 



RDEG = R"5 7 .29 6 
RDDEG=RDOT*5 7 . 296 
DDEG=D*57 . 296 



79 



YAWC=YAWC*57 . 296 
ICOUNT= 1 



C TEST IF WANT TO STOP 
300 IF (TIME . GE . ETIME ) GO TO 400 
C INTEGRATION STEP SIZE BELT 
DELT= 1 . 0 
C INTEGRATION 

U=U+UDOT*DELT 
V = V + VD OT*D ELT 
R = R + RDOT *DELT 
YAW=YAW+R*DELT 
X2=X2+DX2*DELT 

C CONVERT SHIP TO FIXED COORDINATES ON EARTH 
C XDOT = U*DCO S ( YAW ) - V*DS IN ( YAW ) 

C YDOT=U*DSIN (YAW) +V*DCOS ( YAW ) 

C X=X+XDOT*DELT 

C Y=Y+YDOT*DELT 

TIME=TIME+DELT 
I COUNT = I COUNT + 1 
ISE= ISE + LAMD A* YAWE ** 2 
ISR= ISR + D**2 
GO TO 200 

C J= TDIFF = COST FUNCTION 
400 TDIFF = ISE + ISR 

WRITE (6 ,300) ISE , ISR , TDIFF , Kl , TI , T2 
500 FORMAT ( ' f , IX , ’ ISE= ’ , F15 . 7 , ? ISR= ' , F13 . 7 , ' 

IF 13 . 7 ,2X, ' Kl= ? ,F15 . 7 , 2X , 'Tl= ? , F 15 . 7 , 2X , ' T2= ' 
RETURN 
END 

The function minimization subroutine BOXPLX follows 
Then the following two cards must be placed. 

/ / GO . SYSIN DD * 

/* 



TOTALS ' , 
, F15 . 7 ) 
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APPENDIX C 



SYSTEM'S RESPONSE FOR REGULAR SEAS 

//REGURESP JOB (XXXX , XXXX ) , ' RESEARCH ' , CLASS=A 
//-MAIN ORG=NPGVMI.XXXXP 

// EXEC FORTXCG , PARM . FORT= ’ OPT ( 2 ) ' , IMSL=DP , REGION= 1024K 
//FORT. SYS IN DD * 

C IN ORDER TO PERFORM SIMULATION ONLY WHEN GAINS HAVE BEEN 
C OBTAINED CHANGE XS(*) TO X(*) AND DELETE XU(*),AND XL(*) 
COMMON J 
DIMENSION X ( 3 ) 

X(l)=l. 5420017 
X(2)= 141 . 2350922 
X ( 3 ) = 23 . 8 943 6 34 
C CALL PLANT (X) 

C IF ONLY SIMULATION IS WANTED 
CALL PLANT (X) 

WRITE (6,25) 

25 FORMAT (IX,' OPTIMAL GAINS ',/ ) 

DO 30 1=1,3 

30 WRITE(6 ,40)1 ,X(I) 

40 FORMAT (IX, ' X ( ’ ,12, ' )=' , F 14 . 7 ) 

WRITE (6, 50) J 

50 FORMAT ( IX , ’ J = ’,E15.10) 

STOP 

END 

SUBROUTINE PLANT (XX) 

C SUBROUTINE PLANT (XX) SIMULATES THE SHIP 
COMMON TDIFF 
REAL* 8 L,L2 ,L3 ,L4 ,L5 ,L6 

REAL* 8 X , XDOT , Y , YDOT , U , UDOT , V , VDOT , YAW , R , RDOT 
REAL* 8 TIME , ETIME , XUDOT , XUU , XVR , XVV , XDD 
REAL* 8 YV , YR , YD , YVVR , YVRR , YVVV , YRRR , YDDD , YVDOT 
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REAL *8 NV , NR , ND , NVVR , NVRR , NVVV , NRRR , NDDD , NRDOT 
REAL*8 RHO , IZ , FX , FY , MZ , XP , MASS , DELT , MZI , WA , WE 
REAL* 8 DYAWE , YAWE , YAWC , ISE , I SR , LAMDA , D , RYR , RYI , MZR 
REAL* 8 K1 , T1 , T2 , D , X2 , DX2 , S , RX , RY , RZ , TX , TY , TZ , RXR , RXI 
DIMENSION XX ( 3 ) 

C 

C CLOSE LOOP ANALYSIS WITH FILTER 
C 

C INITIAL CONDITIONS FOR INTEGRATION 
C SIMULATION END TIME IN SECONDS 
ETIME= 600 . 

TIME=0 .0 
ICOUNT= 1 

C INITIALIZE THE COST FUNCTION 
ISE=0 . 0 
ISR= 0 . 0 
TDIFF=0 . 0 
LAMDA=8 . 128 

C GAIN COEFFICIENTS TO BE OPTIMIZED 
KI=XX ( 1 ) 

T 1 = XX ( 2 ) 

T 2 = XX ( 3 ) 

C WRITE(6 ,1010) K1,T1,T2 

C1010 FORMAT ( IX , ' K1 =',F15.7,' Tl =',F15.7,’ T2 =',F15.7) 
C X , XDOT , Y , YDOT ARE FIX COORDINATES ON EARTH 
X = 0 . 0 

Y = 0 . 0 
XDOT = 0 . 0 
YDOT=0 . 0 

C U , UDOT , V , VDOT ARE FIX COORDINATES ON SHIP 

V = 0.0 
UDOT=0 . 0 ' 

VDOT=0 .0 
YAW= 0 . 0 
R=0 . 0 
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RDOT=0 . 0 

C ORDERED SPEED IN FEET/ SEC 

C 38.82 FT/ SEC=23 KNOTS 
UC= 38 . 82 

C AT STEADY STATE ACTUAL SPEED (U) = COMMAND SPEED (UC) 

U = UC 

C D = RUDDER ANGLE 
D=0 . 0 
L= 880 . 3 
L2=L**2 
L3=L*L*L 
L4=L*L3 
L5=L*L4 
L6=L*L5 

C SEA DISTURBANCE 

C FORCES IN X,Y DIRECTION COMPUTED IN FORCES 

C MOMENTS IN Z 
FX= 0 . 

FY= 0 . 

MZ = 0 . 

RXR= - . 91037D+03 

RXI=0 .50896D+05 

RYR= - . 2025 6D+ 04 

RYI = . 18077D+06 

MZR= - . 14310D+08 

MZI = - . 16903D+07 

RX = D S QRT ( RXR""" 2 + RX I ** 2 ) 

RY = DSQRT ( RYR ** 2 + R Y I ** 2 ) 

RZ=DSQRT (MZR**2+MZI**2 ) 

TX=DATAN2 (RXI ,RXR) 

TY=DATAN2 (RYI , RYR) 

TZ=DATAN2 (MZI ,MZR) 

C SIGNIFICANT WAVE HEIGHT; SEA STATE 1-0.32,2-0.75,3-2.5, 

•C 4-5.0,5-7.5,6-10.0,7-17.5,8-20.5,9-27.0 
WA= 2 7 . 0 
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C ENCOUNTER FREQUENCY .1,. 2 , .3, .4, .5, .6, .75,1-0, 1.5,2. 5 
WE = 0 . 2 

C HYDRODYNAMIC COEFFICIENTS ARE INSERTED AS PARAMETERS 
RHO=l. 9876 

MASS= ( . 0044)- ( . 5*RHO-L3) 

IZ=(0. 00028 )*( .5*RHO*L5) 

YAWE=0 .0 
X2 = 0 . 0 
DX2 = 0 . 0 
200 CONTINUE 

S = DSQRT(U* * 2 + V * * 2 ) 

C INPUT YAW COMMAND 
YAWC= 0 . 0 

IF (TIME. GE. 0.0) YAWC=0.0 

C ERROR SIGNAL TO DRIVE RUDDER (YAW ACTUAL - YAW ORDERED) 

C ( COMPENSATOR FILTER ) 

YAWE = YAW - YAWC 
DX2= ( YAWE-X2 ) /T2 
ET=K1*(T1*DX2 + X2) 

C AXIAL FORCE HYDRODYNAMIC COEFFICIENTS (SURGE) 

C XUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED FOR 
C DIFFERENT ENCOUNTER ANGLE , SPEED , ENCOUNTER FREQUENCY 
C 

XUDOT = ( - .0001)*( . 5 *RHO*L3 ) 

XU= ( - 0 . 0 25 3 ) * ( . 5 *RHO*L2*S ) 

XUU= ( - 0 . 0 0 0 3 ) * ( . 5 *RHO*L2 ) 

XVR= ( 0 . 0 0 3 9 ) * ( . 5 *RHO*L3 ) 

XVV= ( - . 00 12 )* ( . 5 -RHO-L2 ) 

XDD= ( -0 . 0005 ) * ( . 5*RHO*L2*S**2 ) 

C LATERAL FORCE HYDRODYNAMIC COEFFICIENTS (SWAY) 

C YV = ( - 0 . 00758 )*( . 5 * RHO * L 2 * S ) 

YR= ( 0 . 00 23 ) * ( . 5 *RHO*L3*S ) 

YD= (0 . 00145 )*( . 5 * RHO * L 2 * S * * 2 ) 

YVVR= ( 0 . 0 1 ) * ( . 5 *RHO*L3 / S ) 

YVRR= ( - 0 . 008 ) * ( . 5 - RHO ' ; 'L 4 / S ) 
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YVVV = ( - 0 . 0 3 ) * ( . 5 *RH0*L2 / S ) 

YRRR= ( 0 . 0 0 3 ) * ( . 5 *RH0*L5 / S ) 

YDDD= ( - 0 . 0005) * ( . 5*RHO*L2*S**2 ) 

C YUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED FOR 
C DIFFERENT ENCOUNTER ANGLE , SPEED , ENCOUNTER FREQUENCY 
C 

C YVDOT =(-0.0039) *(.5 *RHO*L3 ) 

C SPEED= 2 3 KNOTS , ENCOUNTER ANGLE= 60 , ENCOUNTER FREQ=0.75 
YVDOT = - . 30908+07 
YV =-0.8127 ID +04 

C MOMENT ABOUT Z-AXIS HYDRODYNAMIC COEFFICIENTS (YAW) 

NV= ( - 0 . 002 13 ) * ( . 5*RHO*L3*S ) 

C NR= ( - 0 . 00 10 5 ) * ( . 5 *RHO*L4*S ) 

ND= ( - 0 . 000 7 ) * ( . 5*RHO*L3*S**2 ) 

NVVR= ( - 0 . 0 15 ) * ( . 5 *RHO*L4 / S ) 

NVRR =(-0.008)*( . 5 ' v RHO“L5 / S ) 

NVVV= ( 0 . 0 1 ) * ( . 5*RHO*L3 / S ) 

NRRR= ( - 0 . 00 6 ) * ( . 5 *RHO*L6 / S ) 

NDDD= ( 0 . 000 1 ) * ( . 5*RHO*L3*S**2 ) 

C NRDOT IS THE ADDED INERTIA TERM WHICH MUST BE CHANGED 
C FOR DIFFERENT ENCOUNTER ANGLE , SPEED , ENCOUNTER FREQUENCY 
C 

C NRDOT = (- 0. 00027) *( .5 *RHO*L5 ) 

C SPEED= 32 KNOTS , ENCOUNTER ANGLE= 60 , ENCOUNTER FREQ=0.75 
NRDOT = -0.2625 ID +12 
NR= - 0 . 53637D+09 
C REGULAR WAVE SEA STATE 

FX=WA*RX*DCOS (WE*TIME + TX) ‘ 

FY= WA“ R Y * DC O S ( WE “ T I ME + T Y ) 

M Z = W A R Z * D C O S ( WE "TIME + TZ ) 

C U ACTUAL SPEED 
C UC COMMANDED SPEED 
C XP = PROPELLER THRUST 
XP=-XUU*UC**2 
C EQUATIONS OF MOTION 
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C UDOT= ( (XVR + MASS)*V*R + XUU*U**2 + XVV*V**2 

C 1 + XDD*D*D + FX. + XP ) / (MASS-XUDOT) 



VDOT= (YV*V + (YR-MASS*U)*R + YD*D + YVVR*V 



2*R 



1 + YVRR*V*R**2 + YVVV*V**3 

2 + YRRR*R**3 + YDDD*D**3 + FY ) / (MASS - YVDOT ) 

RDOT = (NV'-'V + NR*R + ND*D + NVVR*V**2*R + NVRR*V*R**2 
1 + NVW*V**3 + NRRR*R**3 + NDDD*D**3 + MZ ) / (IZ-NRDOT) 
WHEN TO PRINTOUT 

IF (ICOUNT.EQ. 2) GO TO 50 
GO TO 300 

CONVERT RADIANS TO DEGREES 
50 YAWDEG= YAW* 5 7 .296 
RDEG=R*5 7 . 296 



RDDEG=RDOT*57 .296 
DDEG = D*5 7 . 29 6 



YA WC = YA W C * 5 7 .29 6 
WRITE (6,101) TIME , YAWDEG 
101 FORMAT ( IX , F12 . 8 , IX , F12 . 8 ) 
ICOUNT= 1 

C TEST IF WANT TO STOP 
300 IF (TIME . GE . ETIME ) GO TO 400 
C INTEGRATION STEP SIZE DELT 



DELT =1.0 



C INTEGRATION 

U=U+UDOT*DELT 

V = V + VDOT*DELT 
R=R+RDOT*DELT 

Y AW = YAW + R*D ELT 
X2=X2+DX2*DELT 



C CONVERT SHIP TO FIXED COORDINATES ON EARTH 
C XDOT=U*DCOS (YAW ) -V*DSIN(YAW) 

C YDOT=U*DSIN (YAW ) +V*DCOS (YAW ) 

C X = X + XD 0 T *DELT 

C Y=Y+ YDOT*DELT 



TIME = TIME + DELT 
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I COUNT = I COUNT + 1 
ISE=ISE + LAMDA*YAWE**2 
ISR= ISR + D**2 
GO TO 200 

C J = TDIFF = COST FUNCTION 
400 TDIFF = ISE+ISR 

WRITE (6 ,500) ISE , ISR , TDIFF ,K1 , T1 , T2 
500 FORMAT ( ' ’ , IX , ' ISE= ' , F15 . 7 , ’ ISR= ' , F 15 . 7 , ' 

1 F15 . 7 , 2X , ' Kl= ' ,F15.7,2X, ' Tl= ' ,F15.7,2X, ' T2= 
RETURN 
END 

//GO. SYS IN DD * 

/ 



TOTAL = ’ , 
’ , F 15 . 7 ) 
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APPENDIX D 



IRREGULAR SEASTATE FORMULATION 

/ / IRREGAINS JOB (XXXX , XXXX) RESEARCH ' , CLASSIC 
//-MAIN ORG=NPGVMl . XXXXP 

// EXEC FORTXCG , PARM . FORT = ’ OPT(2 ) ' , IMSL=DP , REG I ON = 1024K 
/ / FORT . SYSIN DD * 

C IN ORDER TO PERFORM SIMULATION ONLY WHEN GAINS HAVE BEEN 
C OBTAINED CHANGE XS(*) TO X(*) AND DELETE XU(*) ,AND XL(*) 
DIMENSION X S ( 3 ) , XU ( 3 ) , XL ( 3 ) 

XS(1)=0 . 655751 
XS ( 2 ) = 90 . 548 3 
XS ( 3 ) = 3 6 . 74847 

C XS(I) IS THE STARTING GUESS 

C XL (I) IS THE LOWER LIMIT FOR THE I ’ TH VARIABLE 
C XU (I) IS THE UPPER LIMIT FOR THE I ’ TH VARIABLE 
XL(1)=0.01 
XU ( 1 ) = 2 . 0 
XL (.2) = 20.0 
XU(2) = 180 . 0 
XL ( 3 ) = 5 . 0 
XU ( 3 ) = 180 . 0 

C A DESCRIPTION OF THE FOLLOWING PARAMETERS 
C IS DISCUSSED IN BOXPLX 
R= 9 . / 13 . 

NTA= 1000 
NPR= 100 
NAV= 0 
NV= 3 
IP = 0 

C THE FOLLOWING STATEMENT MUST BE CHANGED TO 
C CALL PLANT (X) 

C IF ONLY SIMULATION IS WANTED 
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CALL BOXPLX ( NV , NAV , NPR , NTA , R , XS , IP , XU , XL , YMN , IER ) 
WRITE (6,25) 

25 FORMAT (IX,’ OPTIMAL GAINS ’,/ ) 

DO 30 1=1,3 

30 WRITE(6 ,40)1 ,XS(I) 

40 FORMAT ( IX , ' X ( ' , 12 , ' ) = ' , F14 . 7 ) 

STOP 

END 

SUBROUTINE PLANT (XX) 

C SUBROUTINE PLANT (XX) SIMULATES THE SHIP 
COMMON TDIFF 
REAL-' 8 L,L2 ,L3 ,L4 ,L5 ,L6 

REAL -8 . X , XDOT , Y , YDOT , U , UDOT , V , VDOT , YAW , R , RDOT 

REAL- 8 TIME , ETIME , XUDOT , XUU , XVR , XVV , XDD 

REAL- 8 YV , YR , YD , YVVR , YVRR , YVVV , YRRR , YDDD , YVDOT 

REAL- 8 NV , NR , ND , NVVR , NVRR , NVVV , NRRR , NDDD , NRDOT 

REAL -8 RHO , IZ , FX , FY , MZ , XP , MASS , DELT 

REAL -8 DYAWE , YAWE , YAWC , ISE , I SR , LAMDA , D 

REAL -8 K1,T1,T2,T3 , T4 , D , X2 , DX2 , X3 ,DX3 ,X4,CH(11) , S 

DIMENSION XX ( 3 ) 

C 

C CLOSE LOOP ANALYSIS WITH FILTER 

C 

C INITIAL CONDITIONS FOR INTEGRATION 

C SIMULATION END TIME IN SECONDS 
ETIME= 600 . 

TIME = 0 . 0 
ICOUNT = 1 

C INITIALIZE THE COST FUNCTION 
ISE = 0 . 0 
ISR= 0 . 0 
TDIFF=0 . 0 
LAMDA= 8 . 128 

C GAIN COEFFICIENTS TO BE OPTIMIZED 
K1=XX(1) 
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T1=XX ( 2 ) 

T2=XX ( 3 ) 

C WRITE(6 , 1010) K1,T1,T2 

C X , XDOT , Y , YDOT ARE FIX COORDINATES ON EARTH 
X=0. 0 
Y=0 . 0 
XDOT= 0 . 0 
YDOT=0 . 0- 

C U , UDOT , V , VDOT ARE FIX COORDINATES ON SHIP 
V=0. 0 
UDOT = 0 . 0 
VDOT = 0 . 0 
YAW= 0 . 0 
R=0 . 0 
RDOT=0 . 0 

C ORDERED SPEED IN FEET/ SEC 

C 38.82 FT/ SEC=23 KNOTS 
UC = 38 . 82 

C AT STEADY STATE ACTUAL SPEED (U) = COMMAND SPEED (UC) 

• U = UC 

C D = RUDDER ANGLE 
D= 0 . 0 
L=880 . 5 
L2=L**2 
L3 =L*L*L 
L4=L*L3 
L5=L*L4 
L6 = L"'L5 

C SEA DISTURBANCE 

C FORCES IN X,Y DIRECTION COMPUTED IN FORCES 

C MOMENTS IN Z 
FX= 0 . 

FY = 0 . 

MZ = 0 . 

C ISEA IS A SWITCH ;ISEA=0 (CALM WATER) ISEA=1 (SEA STATE) 
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ISEA= 1 

C HYDRODYNAMIC COEFFICIENTS ARE INSERTED AS PARAMETERS 
RHO=l. 9876 

MASS = ( . 0044 ) * ( . 5*RHO*L3 ) 

IZ=(0. 00028 )*( . 5*RHO*L5 ) 

YAWE = 0 . 0 
X2 = 0 . 0 
DX2 = 0 . 0 
200 CONTINUE 

S=DSQRT (U**2+V**2 ) 

C INPUT YAW COMMAND 
YAWC=0 . 0 

IF (TIME. GE. 0.0) YAWC=0.0 

C ERROR SIGNAL TO DRIVE RUDDER (YAW ACTUAL - YAW ORDERED) 

C ( CONTROLLER FILTER ) 

YAWE = YAW - YAWC 
DX2= (YAWE-X2 )/T2 
D=K1*(T1*DX2+X2) 

C AXIAL FORCE HYDRODYNAMIC COEFFICIENTS (SURGE) 

C XUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED FOR 
C DIFFERENT ENCOUNTER ANGLE , SPEED , ENCOUNTER FREQUENCY 
C 

XUDOT = ( - . 000 1 ) * ( . 5*RHO*L3 ) 

XU= ( - 0 . 0 2 5 3 ) * ( . 5 *RHO*L2 ■ *S ) 

XUU= (-0.0003) *(.5 *RHO*L2 ) 

XVR= (0.0039) *(.5 *RHO*L3 ) 

XVV= ( - . 00 12 )* ( . 5 "RHO"L2 ) 

XDD= ( - 0 . 00 0 5 ) * ( . 5 *RHO*L2*fc**2 ) 

C LATERAL FORCE HYDRODYNAMIC COEFFICIENTS (SWAY) 

' YV= ( - 0 . 00 7 5 8 ) * ( .5 *RHO*L2 *S ) 

YR- ( 0 . 0 0 2 3 ) * ( . 5 *RHO*L3 * S ) 

YD= ( 0 . 00145 ) * ( . 5*RHO*L2*S**2 ) 

YVVR= ( 0 . 0 1 ) * ( . 5*RHO*L3 / S ) 

YVRR= ( - 0 . 0 0 8 ) * ( . 5 *RHO*L4 / S ) 

YVVV =(-0.03) *(.5 *RHO*L2 / S ) 
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YRRR = ( 0 . 0 0 3 ) * ( . 5 *RHO *L5 / S ) 

YDDD = ( - 0 . 0005 )* ( . 5*RHO*L2*S**2 ) 

C YUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED FOR 
C DIFFERENT ENCOUNTER ANGLE , SPEED , ENCOUNTER FREQUENCY 
C 

C YVDOT = .(-0.0039) *(.5 *RHO*L3 ) 

C SPEED=23 KNOTS, ENCOUNTER FREQUENCY =0.75 
YVDOT = -2304300 . 0 

C MOMENT ABOUT Z-AXIS HYDRODYNAMIC COEFFICIENTS (YAW) 

NV = ( - 0 . 0 0 2 1 3 ) * ( . 5 *RHO*L3 *S ) 

NR= ( - 0 . 00105)* ( . 5*RHO*L4*S ) 

ND = ( - 0 . 0007 ) * ( . 5 * RH 0 * L 3 * S * * 2 ) 

NVVR = ( - 0 . 0 15 ) * ( . 5*RHO*L4 / S ) 

NVRR= ( - 0 . 0 0 8 ) * ( . 5 *RHO*L5 / S ) 

NVVV = ( 0 . 0 1 ) * ( . 5*RHO*L3 / S ) 

NRRR = ( - 0 . 00 6 ) * ( . 5*RHO*L6 / S ) 

NDDD =(0.0001) *(.5 *RHO*L3 * S ** 2 ) 

C NRDOT IS THE ADDED INERTIA TERM WHICH MUST BE CHANGED 
C FOR DIFFERENT ENCOUNTER ANGLE , SPEED , ENCOUNTER FREQUENCY 
C 

C NRDOT= (-0 . 00027 )*( . 5*RHO*L5) 

C SPEED=23 KNOTS, ENCOUNTER FREQUENCY =0.75 
NRDOT = -1.4518E+11 
C SETS SEA STATE TO ZERO 

IF (ISEA. EQ. 1) GO TO 30 
FX= 0 . 

FY=0 . 

MZ = 0 . 

GO TO 35 

C UNIT 12 HAS THE SEA STATE DATA NAMED CH 
C IT MUST BE SYNCHRONIZED BY TIME 
30 READ (12) CH 
FX= CH ( 3 ) 

FY=CH(4 ) 

MZ=CH(8 ) 
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35 CONTINUE 
C U ACTUAL SPEED 
C UC COMMANDED SPEED 
C XP = PROPELLER THRUST 
XP= -XUU*UC**2 



C EQUATIONS OF MOTION 

C UDOT = ( (XVR + MAS S ) *V*R + XUU“U““2 + XW*V**2 

C 1 + XDD*D*D + FX + XP ) / (MASS -XUDOT ) 



VDOT= (YV"V + (YR-MASS*U) *R + YD*D + YVVR*V 



2*R 



C 



C 



1 + YVRR*V*R**2 + YVW*V**3 

2 + YRRR"'R“ " 3 + YDDD*D**3 + 
RDOT= (NV"V + NR*R + ND*D + 

1 + nWV*V**3 + NRRR*R**3 + 
WHEN TO PRINTOUT 

IF (ICOUNT.EQ.il) GO TO 50 
GO TO 300 

CONVERT RADIANS TO DEGREES 
50 YAWDEG= YAW*57.296 



FY )/ (MASS-YVDOT) 
NWR*V**2*R - NVRR*V*R**2 
NDDD*D** 3 + MZ ) / ( I Z - NRDOT ) . 



RDEG=R"57 .296 
RDDEG=RDOT *5 7 . 296 
DDEG=D“57 .296 



YAWC=YAWC*57 .296 
ICOUNT = 1 

C TEST IF WANT TO STOP 
300 IF (TIME . GE . ETIME ) GO TO 400 
C INTEGRATION STEP SIZE DELT 



DELT= 1 . 0 
C INTEGRATION 

U=U+ UDOT “DELT 
V=V + VDOT“DELT 



R=R+RDOT“DELT 
YAW=YAW+R*DELT 
X2=X2+DX2 "DELT 



C CONVERT SHIP TO FIXED COORDINATES ON EARTH 
C XDOT = U"'DCOS (YAW) -V“DSIN (YAW) 
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C YDO T = U*D S I N ( YAW ) + V*DC 0 S ( YAW ) 

C X = X + XDO T “BELT 

C Y = Y + YD 0 T *D ELT 

TIME = TIME +DELT 
I COUNT = I COUNT + 1 
ISE= ISE + LAMD A* Y A WE * * 2 
ISR= ISR + D**2 
GO TO 200 

C J = TDIFF = COST FUNCTION 
400 TDIFF=ISE+ISR 

WRITE (6,500) TDIFF ,K1,T1,T2 
500 FORMAT ( ' IX TOTAL = T ,F15.7,’ K1 =’,F15.7 S 
1 ' T1 =’ ,F15.7,2X, ’T2=’ ,F15.7) 

REWIND 12 

RETURN 

END 

The function minimization subroutine BOXPLX follows 

Then the following three cards must be placed. 

//GO.SYSIN DD * 

/* 

/ GO . FT12F00 1 DD DISP=SHR,DSN=MSS . S2160 . A241 
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APPENDIX E 



SYSTEM'S RESPONSE FOR IRREGULAR SEAS 

//IRRERESP JOB (XXXX , XXXX) RESEARCH CLASS =B 

// -'MAIN ORG-NPGVM1 . XXXXP 

// EXEC FRTXCLGP , IMSL-DP , REGION= 1024K 

//FORT. SYS IN DD * 

C IN ORDER TO PERFORM SIMULATION ONLY WHEN GAINS HAVE BEEN 

C OBTAINED 

DIMENSION XX ( 3 ) 

C OPTIMAL GAINS FOR CONTROLLER 
XX(1)= .9192610 
XX ( 2 ) = 18 .5788116 
XX(3)=9 . 77668 

C THE SUBROUTINE PLANT SIMULATES THE SL-7 CONTAINER SHIP 
CALL PLANT (XX) 

WRITE (6 ,25) 

25 FORMAT (IX, ’OPTIMAL GAINS’ ,/ ) 

DO 30 1=1,3 

30 WRITE (6 ,40)1 ,XX(I) 

40 FORMAT ( IX , ’ XX ( ’ , 12 , ’ ) = ’ , F14 . 7 ) 

STOP 

END 

C 

C SUBROUTINE PLANT (XX) SIMULATES THE SHIP 
SUBROUTINE PLANT (XX) 

COMMON TDIFF 

REAL "8 L,L2 ,L3 ,L4 ,L5 ,L6 

REAL" 8 X , XDOT , Y , YDOT , U , UDOT , V , VDOT , YAW , R , RDOT 
REAL" 8 TIME , ETIME , XUDOT , XUU , XVR , XVV , XDD 
REAL" 8 YV , YR , YD , YVVR , YVRR , YVVV , YRRR , YDDD , YVDOT 
REAL" 8 NV , NR , ND , NVVR , NVRR , NVVV , NRRR , NDDD , NRDOT 
REAL" 8 RHO , IZ , FX , FY , MZ , XP ,MASS , DELT 
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REAL'' 8 DYAWE , YAWE , YAWC , ISE , I SR , LAMDA , D 
REAL “8 K1,T1,T2,D,X2, DX2 , S , CH ( 1 1 ) , DX3 ,X3 , X4 
DIMENSION XX ( 3 ) 

C 

C CLOSE LOOP ANALYSIS WITH FILTER 

C 

INITIAL CONDITIONS FOR INTEGRATION 

C SIMULATION END TIME IN SECONDS 
ETIME= 600 . 

TIME= 0 . 0 
ICOUNT = 1 

C INITIALIZE THE COST FUNCTION 
ISE=0 . 0 
ISR=0 . 0 
TDIFF=0 . 0 
LAMDA = 4 . 2 

C GAIN COEFFICIENTS 
K1=XX ( 1 ) 

T1=XX ( 2 ) 

T 2 = XX ( 3 ) 

C X , XDOT , Y , YDOT ARE FIX COORDINATES ON EARTH 
X=0. 0 
Y= 0 . 0 
XDOT = 0 . 0 
YDOT = 0 . 0 

C U , UDOT , V , VDOT ARE FIX COORDINATES ON SHIP 
V=0.0 
UDOT = 0 . 0 
VDOT = 0 . 0 
YAW=0 . 0 
R= 0 . 0 
RDOT = 0 . 0 
YAW= 0 . 0 

C ORDERED SPEED IN FEET/SEC 

C 38.81 FT / SEC = 23 KNOTS 
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UC= 38 .81 

C AT STEADY STATE ACTUAL SPEED (U) = COMMAND SPEED (UC) 

U = UC 

C D = RUDDER ANGLE 
D=0 . 0 
L = 880 . 5 
L2=L**2 
L3 = L' V L "L 
L4=L*L3 
L5=L*L4 
L6=L*L5 

C SEA DISTURBANCE 

C FORCES IN X,Y DIRECTION . COMPUTED IN FORCES 
C MOMENTS IN Z 
FX = 0 . 

FY = 0 . 

MZ = 0 . 

C I SEA IS A SWITCH; ISEA=0 (CALM WATER) ISEA=1 (SEA STATE) 
ISEA= 1 

C HYDRODYNAMIC COEFFICIENTS ARE INSERTED AS PARAMETERS 
RHO=l. 9876 

MAS S = ( . 0044 ) * ( . 5*RHO*L3 ) 

I Z = ( 0 . 0 0 0 2 8 ) * ( . 5 *RHO*L5 ) 

YAWE = 0 . 0 
X2=0 .0 
DX2=0 . 0 
X3 = 0 . 0 
DX3 = 0 . 0 
X4 = 0 . 0 

200 CONTINUE 

S=DSQRT(U**2 + V**2) 

C INPUT YAW COMMAND 
YAWC = 0 . 0 

IF (TIME. GE. 0.0) YAWC=0.0 

C ERROR SIGNAL TO DRIVE RUDDER (YAW ACTUAL - YAW ORDERED) 
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C ( COMPENSATOR FILTER ) 

YAWE= YAW - YAWC 
DX2= (YAWE-X2 ) /T2 
X4 = K 1 * ( T 1*DX2 + X2 ) 

DX3= (X4-X3 )/T4 
D= (T3*DX3+X3) 

C AXIAL FORCE HYDRODYNAMIC COEFFICIENTS (SURGE) 

C 

C XUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED FOR 
C DIFFERENT ENCOUNTER ANGLE AND SPEED. 

C XUDOT = ( - . 0 0 0 1 ) * ( . 5 *RHO*L 3 ) 

XUU= ( - 0 . 0 0 0 3 ) * ( . 5 *RHO*L2 ) 

XVR= ( 0 . 00 3 9 ) * ( . 5*RHO*L3 ) 

XVV= ( - . 00 12 ) * ( . 5*RHO*L2 ) 

XDD= ( -0 . 0005)* ( . 5*RHO*L2*S**2 ) 

C LATERAL FORCE HYDRODYNAMIC COEFFICIENTS (SWAY) 

YV = ( - 0 . 0 0 7 5 8 ) * ( .5 *RHO*L2 * S ) 

YR= (0 . 0023 ) * ( . 5*RHO*L3*S ) 

YD= (0 . 00145 )* ( . 5*RHO*L2*S**2 ) 

YVVR= ( 0 . 0 1 ) * ( . 5 *RHO*L3 / S ) 

YVRR= ( - 0 . 0 0 8 ) * ( . 5 *RH0*L4 / S ) 

YWV = ( - 0 . 0 3 ) * ( . 5 *RHO*L2 / S ) 

YRRR= ( 0 . 0 0 3 ) * ( . 5 *RHO*L5 / S ) 

YDDD= ( - 0 . 0005 ) * ( . 5*RHO*L2*S**2 ) 

C YUDOT IS THE ADDED MASS TERM WHICH MUST BE CHANGED FOR 
C DIFFERENT ENCOUNTER ANGLE AND SPEED. 

C 

C YVDOT = ( - 0 . 00 3 9 ) * ( . 5 *RHO*L3 ) 

YVDOT= -3654800 . 00 

C MOMENT ABOUT Z-AXIS HYDRODYNAMIC COEFFICIENTS (YAW) 

NV= ( - 0 . 00 2 1 3 ) * ( . 5 *RHO*L3 *S ) 

NR= ( - 0 . 00 105 )* ( . 5*RHO*L4*S ) 

ND= ( - 0 . 0007 ) * ( . 5 *RHO *L 3 * S ** 2 ) 

NVVR= ( - 0 . 0 15 ) * ( . 5 *RHO*L4 / S ) 

NVRR= ( - 0 . 008 ) * ( . 5*RHO*L5 / S ) 
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NVVV = ( 0 . 0 1 ) * ( . 5 *RH0*L3 / S ) 

NRRR = ( - 0 . 00 6 ) * ( . 5 *RHO*L6 / S ) 

NDDD= (0 . 000 1 )* ( . 5 * RHO * L 3 * S * * 2 ) 

C NRDOT IS THE ADDED INERTIA TERM WHICH MUST BE CHANGED 
C FOR DIFFERENT ENCOUNTER ANGLE AND SPEED. 

C 

C NRDOT = (-0 . 00027 )*( . 5*RHO*L5 ) 

NRDOT = - 2 . 1815E+11 
C SETS SEA STATE TO ZERO 

IF (ISEA.EQ. 1) GO TO 30 
FX= 0 . 



FY = 0 . 
MZ = 0 . 



GO TO 35 

C UNIT 12 HAS THE SEA STATE DATA NAMED CH 
C IT MUST BE SYNCHRONIZED BY TIME 
30 READ (12) CH 

FX= CH ( 3 ) 

FY= CH (4 ) 

MZ= CH (8) 

35 CONTINUE 

C U ACTUAL SPEED 
C UC COMMANDED SPEED 
C XP = PROPELLER THRUST 
X P = - XUU * U C * * 2 
C EQUATIONS OF MOTION 

C UDOT = ( (XVR + MASS ) *V*R + XUU*U**2 + XVV*V**2 

C 1 + XDD*D*D + FX + XP ) / (MASS - XUDOT ) 

VDOT = (YV*V + ( YR - MAS S * S ) *R + YD*D + YWR*V**2*R 



C 



1 + Y V R R * V * R * * 2 + YVW*V**3 
1 + YRRR*R**3 + YDDD*D**3 + 
RDOT= (NV*V + NR*R + ND*D + 
1 + NVVV* V** 3 + NRRR*R**3 + 
WHEN TO PRINTOUT 

IF (ICOUNT.EQ.2 ) GO TO 50 



FY ) / (MASS- YVDOT) 
NWR*V**2*R + NVRR*V*R**2 
NDDD*D**3 + MZ ) / ( I Z - NRDOT ) 
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GO TO 300 

C CONVERT RADIANS TO DEGREES 
5 0 YAWDEG = YAW*5 7.296 
RDEG=R*57 .296 
RDDEG=RDOT*57 . 296 
DDEG=D"57 .296 
YAWC=YAWC*57 . 296 
WRITE (6,100) TIME, YAWDEG 
100 FORMAT ( IX , F 12 . 8 , IX , F 12 . 8 ) 

ICOUNT = 1 

C TEST IF WANT TO STOP 
300 IF (TIME . GE . ETIME ) GO TO 400 

C INTEGRATION STEP SIZE DELT 
DELT= 1 . 

C INTEGRATION 

U=U+UDOT*DELT 
V = V + VD 0 T ' V D E LT 
R=R + RDOT' v DELT 
YAW = YAW + R' v DELT 
X2=X2+DX2*DELT 
X3 = X3 +DX3' V DELT 

C CONVERT SHIP TO FIXED COORDINATES ON EARTH 
XDOT =U*DCO S ( YAW ) - V* DS IN ( YAW ) 

YDOT=U*DSI N ( YAW ) + V*DCO S ( YAW ) 

X= X + XDOT' v DELT 
Y=Y+YDOT*DELT 
TIME=TIME+DELT 
I COUNT = I COUNT + 1 
ISE=ISE + LAMDA*YAWE**2 
ISR= ISR + D**2 
GO TO 200 

C J = TDIFF = COST FUNCTION 
400 TDIFF = ISE+ISR 

WRITE (6 , 500 ) ISE, ISR, TDIFF 
500 FORMAT ( ' 1 ’ ,5X, ' ISE= ' , F15 . 7 , ' ISR=',F15.7, 



100 



1 ’ TOTAL= ' ,F15 . 7) 

STOP 

END 

//GO. SYS IN DD * 

/* 

/ /GO . FT12F00 1 DD DISP = SHR , DSN = MSS . S2 160 . A2 11 
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APPENDIX F 



MODIFIED MINIMIZATION SUBROUTINE 

SUBROUTINE BOXPLX (NV , NAV , NPR , NTZ , RZ , XS , IP , BU , BL , YMN , IER ) 
C 

DIMENSION V(50 ,50) , FUN(50), SUM(25), CEN(25), XS(NV), 
1BU(NV) , BL ( NV ) 

C 

KV = 5 
EP = l.E-4 
NTA = 2000 

IF (NTZ.GT.O) NTA = NTZ 
R = RZ 

IF (R.LE.O. .OR.R.GE. 1. ) R=l./3. 

NVT = NV+NAV 
C 

C TOTAL VARS, EXPLICIT PLUS IMPLICIT 

NT = 0 

C CURRENT TRIAL NO. 

NPT = 0 

C CURRENT NO. OF PERMISSIBLE TRIALS 

NTFS = 0 

C CURRENT NO. OF TIMES F HAS BEEN ALMOST UNCHANGED 
C 

C CHECK FEASIBILITY OF START POINT 

C 

DO 4 1=1, NV 
VT = XS(I) 

IF ( BL ( I ) .LE. VT) GO TO 1 
II = -I 
VT = BL ( I ) 

GO TO 2 

1 IF ( BU ( I ) .GE. VT) GO TO 3 



102 



II = I 
VT = BU (I ) 

2 IF (NPR.GT.O) WRITE (6,49) II 

3 V(I,1) = VT 
CEN(I) = VT 

IF (IP.EQ. 1) GO TO 4 

BL ( I ) = BL ( I ) + AMAXl ( EP , EP*ABS ( BL ( I ) ) ) 

BU ( I ) = BU(I)-AMAX1(EP , EP-'ABS (BU ( I ) ) ) 

4 SUM (I) = VT 
C 

c 

NCE = 1 

C NUMBER OF CONSTRAINT EVALUATIONS 
1 = 1 

IF (KE(V(1, 1) ) .EQ.O) GO TO 5 
IF (NPR.LE.O) GO TO 12 
WRITE (6,50) 

GO TO 12 

5 NFE = 1 
C 

C NUMBER OF VERTICES (K) = 2 TIMES NO. OF VARIABLES. 

K = (2*NV)/3 
C 

C NUMBER OF DISPLACEMENTS ALLOWED. 

NLIM = 5 *NV +10 
C 

C NUMBER OF CONSECUTIVE TRIALS WITH UNCHANGED FE TO 
C TO TERMINATE 

NCT = NLIM+NV 
ALPHA =1.3 
FK = K 
FKM = FK-1. 

BETA = ALPHA+1. 

C 

C INSURE SEED OF RANDOM NUMBER GENERATOR IS ODD. 
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IQR = R* 1 . E 7 

IF (MOD (IQR, 2) . EQ. 0) IQR=IQR+101 
C 

C SET UP INITIAL VERTICES 
FUN ( 1 ) = FE ( V ( 1 , 1 ) ) 

YMN = FUN ( 1 ) 

6 FI = 1. 

FUNOLD = FUN ( 1 ) 

C 

DO 15 1=2, K 
FI = FI+1. 

LIMT = 0 

7 LIMT = LIMT + 1 
C 

C END CALCULATION IF FEASIBLE CENTROID CANNOT BE FOUND. 

IF (LIMT.GE.NLIM) GO TO 11 
C 

DO 8 J= 1 , NV 
C 

C RANDOM NUMBER GENERATOR (RANDU) 

IQR = IQR- 6553 9 

IF (IQR.LT.O) IQR = IQR+2147483647+ 1 
RQX = IQR 

RQX = RQX- .4656613E-9 

V(J,I) = BL(J)+RQX*(BU( J)-BL(J) ) 

IF (IP.EQ.l) V(J,I)=AINT(V(J,I)+ .5) 

8 CONTINUE 
C 

DO 10 L= 1 , NLIM 
NCE = NCE+ 1 

IF (KE ( V ( 1 , I ) ) . EQ . 0 ) GO TO 13 
C 

DO 9 J = 1 , NV 

VT = . 5 - ( V ( J , I ) + CEN ( J ) ) 

IF (IP.EQ.l) VT = AINT (VT+ . 5 ) 
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V(J,I) = VT 
9 CONTINUE 
C 

10 CONTINUE 
C 

11 IF (NPR.LE.O) GO TO 12 
WRITE (6,51) I 

CALL BOUT (NT , NPT , NFE , NCE , NV , NVT , V , I , FUN , CEN , I ) 

12 IER = -1 
GO TO 48 

C 

13 DO 14 J = 1 , NV 

SUM ( J ) = SUM (J)+V(J,I) 

14 CEN(J) = SUM(J)/FI 
C 

C TRY TO ASSURE FEASIBLE CENTROID FOR STARTING. 

NCE = NCE+1 

IF (KE(CEN) .EQ.O) GO TO 60 
SUM ( J ) = SUM ( J ) -V(J,I) 

GO TO 7 

60 NFE = NFE+ 1 

FUN ( I ) = FE ( V ( 1 , I ) ) 

15 CONTINUE 
C 

C END OF LOOP SETTING OF INITIAL COMPLEX. 

IF (NPR.LE.O) GO TO 17 

CALL BOUT (NT, NPT, NFE, NCE, NV, NVT, V,K, FUN, CEN, 0) 
C 

C FIND THE WORST VERTEX, THE ' J ’ TH . 

J = 1 
C 

DO 16 1 = 2, K 

IF ( FUN ( J ) . GE . FUN ( I ) ) GO TO 16 
J = I 

16 CONTINUE 
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c 

C BASIC LOOP. ELIMINATE EACH WORST VERTEX IN TURN. 

C IT MUST BECOME NO LONGER WORST, NOT MERELY IMPROVED. 
C FIND NEXT -TO -WORST VERTEX, THE ' JN ' TH ONE. 

17 JN = 1 

IF (J.EQ. 1) JN = 2 
C 

DO 18 1=1, K 
IF (I.EQ.J) GO TO 18 
IF (FUN(JN) .GE.FUN(I) ) GO TO 18 
JN = I 

18 CONTINUE 
C 

C LIMT = NUMBER OF MOVES DURING THIS TRIAL TOWARD THE 
C CENTROID DUE TO FUNCTION VALUE. 

LIMT = 1 
C 

C COMPUTE CENTROID AND OVER REFLECT WORST VERTEX. 

C 

DO 19 1=1, NV 
VT = V ( I , J ) 

SUM ( I ) = SUM ( I ) - VT 

CEN(I) = SUM ( I ) / FKM 

VT = BETA“CEN ( I ) - ALPHA* VT 

IF (IP.EQ.l) VT = AINT (VT+ . 5 ) 

C 

C INSURE THE EXPLICIT CONSTRAINTS ARE OBSERVED. 

19 V(I,J) = AMAX1 ( AMINl ( VT , BU ( I ) ) , BL ( I ) ) 

C 

NT = NT + 1 
C 

C CHECK FOR IMPLICIT CONSTRAINT VIOLATION. 

C 

20 DO 25 N= 1 , NLIM 
NCE = NCE+1 



106 



IF (KE ( V ( 1 , J ) ) . EQ . 0 ) GO TO 26 
C 

C EVERY 'KV'TH TIME, OVER- REFLECT THE OFFENDING VERTEX 
C THROUGH THE BEST VERTEX. 

IF (MOD (N , KV ) . NE . 0 ) GO TO 22 
CALL FBV (K , FUN ,M) 

C 

DO 21 1=1, NV 

VT = BETA*V ( I , M ) - ALPHA “V ( I , J ) 

IF (IP.EQ.l) VT = AINT (VT+ . 5 ) 

21 V(I,J) = AMAX 1 ( AMIN 1 ( VT , BU ( I ) ) , BL ( I ) ) 

C 

GO TO 24 
C 

C CONSTRAINT VIOLATION: MOVE NEW POINT TOWARD CENTROID. 

C 

22 DO 23 1=1, NV 

VT = . 5*(CEN(I)+V(I , J) ) 

IF (IP.EQ.l) VT = AIN.T (VT+ . 5 ) 

V (I , J ) = VT 

23 CONTINUE 
C 

24 NT = NT + 1 

25 CONTINUE 
C 

IER = 1 
C 

C CANNOT GET FEASIBLE VERTEX BY MOVING TOWARD CENTROID, 

C OR BY OVER-REFLECTING THRU THE BEST VERTEX. 

IF (NPR.LE.O) GO TO 42- 
WRITE (6,52) NT , J 

CALL BOUT ( NT , NPT , NFE , NCE , NV , NVT , V , K , FUN , CEN , J ) 

GO TO 42 
C 

C FEASIBLE VERTEX FOUND, EVALUATE THE OBJECTIVE FUNCTION. 
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26 NFE = NFE+1 
FUNTRY = FE (V ( 1 , J ) ) 

C 

C TEST TO -SEE IF FUNCTION VALUE HAS NOT CHANGED. 

AFO = ABS( FUNTRY- FUNOLD) 

AMX = AMAX1 (ABS(EP -'FUNOLD ) ,EP) 

C 

C ACTIVATE THE FOLLOWING TWO STATEMENTS 
C FOR DIAGNOSTIC PURPOSES ONLY. 

C WRITE (6,99) J , AFO , AMX , FUNTRY , FUNOLD , FUN (J ), FUN ( JN ) , 

C 1NTFS , N 

C 99 FORMAT ( IX , 13 , 6E15 . 7 , 215 ) 

IF (AFO. GT. AMX) GO TO 27 
NTFS = NTFS + 1 
IF (NTFS .LT.NCT) GO TO 28 
IER = 0 

IF (NPR.LE.O) GO TO 42 
WRITE (6,53) K 

CALL BOUT ( NT , NPT , NFE , NCE , NV , NVT , V , K , FUN , CEN , 0 ) 

GO TO 42 

27 NTFS = 0 
C 

C IS THE NEW VERTEX NO LONGER WORST? 

28 IF ( FUNTRY.LT. FUN ( JN ) ) GO TO 34 
C 

C TRIAL VERTEX IS STILL WORST; ADJUST TOWARD CENTROID. 

C EVERY ’KV'TH TIME, OVER-REFLECT THE OFFENDING VERTEX 

C THROUGH THE BEST VERTEX. 

LIMT = LIMT + 1 

IF (MOD (LIMT, K.V ) . NE . 0 ) GO TO 30 
CALL FBV (K , FUN , M) 

C 

DO 29 1=1, NV 

VT = BETA*V ( I , M ) - ALP HA* V ( I , J ) 

IF (IP.EQ.l) VT = AINT ( VT + . 5 ) 
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29 V(I,J).= AMAX1 (AMIN1 (VT , BU (I)),BL(I)) 

C 

GO TO 32 
C 

30 DO 31 I = 1 , NV 

VT = . 5*(CEN(I)+V(I , J) ) 

IF (IP.EQ.l) VT = AINT ( VT + . 5 ) 

V(I,J) = VT 

31 CONTINUE 
C 

32 IF (LIMT.LT.NLIM) GO TO 33 
C 

C CANNOT MAKE THE ' J ' TH VERTEX NO LONGER WORST BY 
C DISPLACING TOWARD THE CENTROID OR BY OVER- REFLECTING 
C THRU THE BEST VERTEX. 

IER = 2 

IF (NPR . LE. 0) GO TO 42 
WRITE (6,52) NT, J 

CALL BOUT ( NT , NPT , NFE , NCE , NV , NVT , V , K , FUN , CEN , J ) 

GO TO 42 

33 NT = NT + 1 • 

GO TO 20 

C 

C SUCCESS: WE HAVE A REPLACEMENT FOR VERTEX J. 

34 FUN ( J ) = FUNTRY 
FUNOLD = FUNTRY 
NPT = NPT + 1 

C 

C STOP AT THE 100 'TH PERMISSIBLE TRIAL 
IF (MOD(NPT , 100) .EQ. 0) GO TO 48 
C 

DO 36 1=1, NV 
SUM ( I ) = 0. 

c 

DO 35 N=1,K 
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35 SUM ( I ) = SUM (I)+V(I,N) 

C ■ 

CEN(I) = SUM ( I ) / FK 

36 CONTINUE 
C 

LC = 0 
GO TO 39 
C 

37 DO 38 1=1, NV 

38 SUM ( I ) = SUM (I)+V(I,J) 

C 

LC = J 
C 

39 IF (NPR.LE.O) GO TO 40 

IF (MOD (NPT , NPR) . NE . 0 ) GO TO 40 
C 

CALL BOUT ( NT , NPT , NFE , NCE , NV , NVT , V , K , FUN , CEN , LC ) 

C 

C HAS THE MAX. NUMBER OF TRIALS BEEN REACHED WITHOUT 
C CONVERGENCE ? IF NOT, GO TO NEW TRIAL. 

40 IF (NT.GE.NTA) GO TO 41 
C 

C NEXT-TO-WORST VERTEX NOW BECOMES WORST. 

J = JN 
GO TO 17 

41 IER = 3 

IF (NPR.GT.O) WRITE (6,54) 

C 

C COLLECTOR POINT FOR ALL ENDINGS. 



c 


1) 


CANNOT 


DEVELOP 


FEASIBLE VERTEX. 




IER = 1 


c 


2) 


CANNOT 


DEVELOP 


A NO-LONGER- WORST VERTEX. 


IER = 2 


c 


3) 


FUNCTION VALUE 


UNCHANGED FOR K 


TRIALS . 


IER = 0 


c 


4) 


LIMIT 


ON TRIALS 


REACHED . 




IER = 3 


c 


5) 


CANNOT 


FIND FEASIBLE VERTEX AT 


START . 


IER =-l 



42 CONTINUE 
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c 

C FIND BEST VERTEX. 

CALL FBV ( K , FUN , M ) 

IF (IER.GE.3) GO TO 44 

C RESTART IF THIS SOLUTION IS SIGNIFICANTLY BETTER THAN 
C THE PREVIOUS OR IF THIS IS THE FIRST TRY. 

IF (NPR.LE.O) GO TO 43 
WRITE (6,55) (M , YMN , FUN (M) ) 

43 IF (FUN(M) .GE. YMN) GO TO 47 

IF (ABS(FUN(M)-YMN) . LE . AMAX1 ( EP , EP*YMN ) ) GO TO 47 
C GIVE IT ANOTHER TRY UNLESS LIMIT ON TRIALS REACHED. 

44 YMN = FUN (M) 

FUN ( 1 ) = FUN (M) 

C 

DO 45 1=1, NV 
CEN(I) = V ( I , M ) 

SUM ( I ) = V ( I , M ) 

45 V(I,1) = V ( I , M ) 

C 

DO 46 1=1, NVT 

46 XS(I) = V ( I ,M) 

C 

IF (IER.LT.3) GO TO 6 

47 IF (NPR.LE.O) GO TO 48 

CALL BOUT (NT , NPT , NFE , NCE , NV , NVT , V ,K , FUN , V ( 1 ,M) , - 1 ) 
WRITE (6,56) FUN (M) 

48 RETURN 
C 

49 FORMAT ( 50H0INDEX AND DIRECTION OF OUTLYING 
1VARIABLE AT STARTI5 ) 

50 FORMAT ( 5 OHOIMPLICIT CONSTRAINT VIOLATED AT START. 
1DEAD END . ) 

51 FORMAT (’OCANNOT FIND FEASIBLE ’, I 4 ,’ TH VERTEX OR 
1CENTROID AT START . ’ ) 

52 FORMAT ( 10H0AT TRIAL I4,54H CANNOT FIND FEASIBLE 
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1VERTEX WHICH IS NO 

2LONGER WORST, 14, 15X, ’RESTART FROM BEST VERTEX.') 

53 FORMAT (40H0FUNCTION HAS BEEN ALMOST UNCHANGED 
IFOR 15 , 7H TRIALS) 

54 FORMAT (27H0LIMIT ON TRIALS EXCEEDED. ) 

55 FORMAT ( ’ OBEST VERTEX IS NO. ’,13,’ OLD MIN WAS ’, 

1E15.7, ' NEW MIN IS ’,E15.7) 

56 FORMAT ( ’ OMIN OBJECTIVE FUNCTION IS ' ,E15.7) 

END 

SUBROUTINE FBV (K , FUN ,M) 

DIMENSION FUN (50) 

M = 1 
C 

DO 1 1 = 2, K 

IF (FUN(M) .LE.FUN(I) ) GO TO 1 
M = I 

1 CONTINUE 
C 

RETURN 

END 

SUBROUTINE BOUT (NT , NPT , NFE , NCE , NV , NVT , V , K , FN , C , IK ) 
DIMENSION V(50 ,50) , FN(50), C(25) 

WRITE (6,4) NT, NPT, NFE, NCE 
C 

DO 1 1=1, K 

WRITE (6,5) FN(I) , (V(J, I) , J=1 ,NV) 

IF (NVT.LE.NV) GO TO 1 
NVP = NV + 1 

WRITE (6,6) (V(J,I) , J=NVP,NVT) 

1 CONTINUE 
C 

IF (IK.NE.O) GO TO 2 
C 

WRITE (6,7) (C(I) ,I=1,NV) 

RETURN 
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2 IF (IK.GE.O) GO TO 3 
WRITE (6,8) (C(I) ,I=1,NV) 

RETURN 

3 WRITE (6,9) IK, ( C ( I ) , I = 1 , NV ) 

RETURN 

4 FORMAT ('ONO. TOTAL TRIALS = ' ,I5,4X,'NO. FEASIBLE 
1TRIALS = ' ,' 

215 , 4'X , ' NO . FUNCTION EVALUATIONS = ’ ,I5,4X,'NO. 
3CONSTRAINT EVALUATIONS 

4= ' , 15 / ' 0 FUNCTION VALUE' , 6X, ' INDEPENDENT 

5VARIABLES/ DEPENDENT 

60R IMPLICIT CONSTRAINTS') 

5 FORMAT (1H , EI8 . 7 , 2X , 7E14 . 7 / (21X , 7E14 . 7 ) ) 

6 FORMAT ( 2 IX , 7EI4 . 7 ) 

7 FORMAT ( 10H0CENTROID I IX , 7E14 . 7/ (21X, 7EI4 . 7 ) ) 

8 FORMAT (’0 BEST VERTEX ' , 7X , 7E14 . 7 / ( 2 IX , 7E 14 . 7 ) ) 

9 FORMAT ( ' OCENTROID LESS VX’ ,12 ,2X, 7EI4. 7/ (2IX, 7EI4 . 7 ) ) 
END 

FUNCTION FE (X) 

DIMENSION X ( 3 )• 

COMMON TDIFF 
CALL PLANT (X) 

FE=TDIFF 

RETURN 

END 

FUNCTION KE (X) 

DIMENSION X ( 3 ) 

KE = 0 

RETURN 

END 
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